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Abstract 

Heat  stroke  (HS)  is  a  life-threatening  illness  induced  by  prolonged  exposure  to  a  hot  environment  that  causes  central 
nervous  system  abnormalities  and  severe  hyperthermia.  Current  data  suggest  that  the  pathophysiological  responses  to  heat 
stroke  may  not  only  be  due  to  the  immediate  effects  of  heat  exposure  per  se  but  also  the  result  of  a  systemic  inflammatory 
response  syndrome  (SIRS).  The  observation  that  pro-  (e.g.,  IL-1)  and  anti-inflammatory  [e.g.,  IL-10)  cytokines  are  elevated 
concomitantly  during  recovery  suggests  a  complex  network  of  interactions  involved  in  the  manifestation  of  heat-induced 
SIRS.  In  this  study,  we  measured  a  set  of  circulating  cytokine/soluble  cytokine  receptor  proteins  and  liver  cytokine  and 
receptor  mRNA  accumulation  in  wild-type  and  tumor  necrosis  factor  (TNF)  receptor  knockout  mice  to  assess  the  effect  of 
neutralization  of  TNF  signaling  on  the  SIRS  following  HS.  Using  a  systems  approach,  we  developed  a  computational  model 
describing  dynamic  changes  (intra-  and  extracellular  events)  in  the  cytokine  signaling  pathways  in  response  to  HS  that  was 
fitted  to  novel  genomic  (liver  mRNA  accumulation)  and  proteomic  (circulating  cytokines  and  receptors)  data  using  global 
optimization.  The  model  allows  integration  of  relevant  biological  knowledge  and  formulation  of  new  hypotheses  regarding 
the  molecular  mechanisms  behind  the  complex  etiology  of  HS  that  may  serve  as  future  therapeutic  targets.  Moreover,  using 
our  unique  modeling  framework,  we  explored  cytokine  signaling  pathways  with  three  in  silico  experiments  (e.g.  by 
simulating  different  heat  insult  scenarios  and  responses  in  cytokine  knockout  strains  in  silico). 
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Introduction 

Heat  stroke  (HS)  is  a  life-threatening  illness  characterized  by 
profound  central  nervous  system  dysfunction,  severely  elevated 
core  temperature,  as  well  as  organ  and  tissue  damage  resulting 
from  environmental  heat  exposure  [1].  Environmental  heat 
exposure  is  one  of  the  most  deadly  natural  hazards  in  the  United 
States  with  ~200  deaths  per  year.  In  the  past  two  decades, 
extreme  heat  exposure  claimed  more  American  lives  than  the 
combined  effects  of  hurricanes,  lightning,  earthquakes,  floods  and 
tornadoes  [2].  HS  is  also  an  international  hazard  as  demonstrated 
by  the  high  incidence  of  death  (>15,000  individuals)  during  the 
2003  heat  wave  in  France  [3,4],  Clinical  and  experimental 
evidence  suggests  that  the  pathophysiological  responses  to  HS 
are  the  result  of  a  systemic  inflammatory  response  syndrome 
(SIRS)  that  ensues  following  HS  collapse.  The  SIRS  is  regarded  as 
a  response  to  bacteria  and/ or  endotoxin  leakage  across  ischemic- 
damaged  gut  epithelial  barrier  membranes,  which  stimulates 
cytokine  and  other  inflammatory  pathways  that  are  thought  to 
mediate  a  variety  of  pathophysiological  responses.  The  liver  has 
been  implicated  as  an  early  key  player  in  the  heat-induced  SIRS 


based  on  its  function  as  a  major  site  of  endotoxin  clearance  [5]. 
Cytokines  are  important  regulators  of  the  acute-phase  response 
(APR)  to  inflammation/injury  and  have  been  implicated  as 
mediators  of  the  SIRS  with  HS  [6], 

Several  studies  have  characterized  peripheral  cytokine  distur¬ 
bances  in  HS  patients.  At  the  time  of  clinical  admission  or  shortly 
after  cooling,  the  concentration  of  circulating  interleukin  (IL)-la, 
IL-1 /),  IL-1  receptor  antagonist  (IL-IRa),  IL-6,  soluble  IL-6 
receptor  (sIL-6R),  IL-10,  interferon  (IFN)y,  tumor  necrosis  factor 
(TNFja,  and/or  soluble  TNF  receptors  subtype  I  (sTNF-RI)  and 
subtype  II  (sTNF-RII)  have  been  shown  to  be  elevated  in  some  HS 
patients  [7-10].  Unfortunately,  circulating  cytokines  are  often 
determined  primarily  at  end-stage  HS,  which  has  limited  our 
understanding  of  the  time  course  of  changes  in  the  balance  of  these 
mediators  during  progression  of  the  SIRS.  Moreover,  the  complex 
interactions  among  cytokines  that  mediate  the  APR  and  SIRS 
remain  unknown.  Development  of  a  conscious  mouse  model  that 
simulates  the  human  pathophysiological  responses  to  HS  has 
demonstrated  that  plasma  concentrations  of  IL- 1  /?,  IL-6,  IL-10, 
and  IL-12p40  are  increased  in  a  time-  and  core  temperature  (In¬ 
dependent  manner  [6].  Furthermore,  concomitant  elevation  of 
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pro-  (e.g.,  IL-1)  and  anti-inflammatory  (e.g.,  IL-10)  cytokines 
suggests  that  a  complex  network  of  interactions  orchestrates  the 
SIRS  during  HS  recovery. 

Accompanying  elevations  in  cytokines,  organ  (kidney,  liver, 
spleen)  and  tissue  (gut,  skeletal  muscle)  damage  are  common 
manifestations  of  the  HS  syndrome  [11—13].  The  liver  is  a  major 
immune  organ  known  to  produce  and  respond  to  cytokines  during 
inflammation  [14]  and  damage  to  this  organ  is  primarily  observed 
in  long-term  survivors  of  HS  [15].  However,  it  is  unknown  if  liver 
damage  is  a  consequence  of  direct  thermal  injury  or  cytokine- 
induced  pathophysiological  changes  associated  with  the  SIRS, 
indicating  the  importance  of  correlating  changes  in  circulating 
cytokine  levels  with  inflammatory  changes  occurring  at  the  organ 
and/or  tissue  level.  Helwig  and  Leon  [14]  determined  plasma, 
liver,  and  spleen  mRNA  accumulation  patterns  for  the  IL- 1  family 
members  in  mice  following  HS;  increased  IL- la,  IL-1  [i,  and  IL- 1 
receptor  subtype  I  (IL-1RI)  and  subtype  II  (IL-1RII)  mRNA 
accumulation  in  the  liver  and  spleen  suggested  these  organs  may 
contribute  to  circulating  IL-1  family  protein  levels  following  HS, 
but  the  absence  of  studies  on  protein  translation  that  include 
protein  tagging  precluded  a  conclusive  association. 

TNF-a  has  been  shown  to  have  deleterious  actions  in 
endotoxemia  [16,17]  and  it  has  been  assumed  that  this  cytokine 
functions  similarly  in  HS.  Leon  [18]  conducted  the  first  studies 
using  TNF  receptor  subtype  I  (TNF-RI)  and  subtype  II  (TNF-RII) 
knockout  mice  (TNFR  KO;  i.e.,  can  produce  TNF,  but  do  not 
have  signaling  receptors)  and  showed  slower  heating  and  faster 
cooling  rates  in  KO  compared  to  wild-type  mice.  Although  the  Tc 
responses  displayed  by  TNFR  KO  mice  would  be  considered 
protective  against  HS,  these  mice  showed  a  trend  towards 
increased  mortality  compared  to  their  wild-type  controls  during 
the  second  day  of  recovery  (40%  vs.  100%  survival,  respectively). 
This  preliminary  study  indicates  that  TNF  might  have  time- 
dependent  pro-  (early)  and  anti-inflammatory  (late)  actions  in  the 
HS  syndrome,  although  the  mechanisms  mediating  the  early 
actions  of  this  cytokine  in  HS  remain  unidentified.  It  is  important 
to  elucidate  the  pro-inflammatory  actions  of  TNF  in  the  heat- 
induced  SIRS  to  determine  if  this  cytokine  may  be  an  important 
therapeutic  target  to  mitigate  morbidity/ mortality  associated  with 
this  syndrome. 

Using  a  conscious  mouse  HS  model,  we  showed  previously  that 
several  pro-  and  anti-inflammatory  cytokines  are  elevated  in  the 
circulation  and  liver  following  HS  collapse  [14,18].  TNF  is  known 
to  interact  with  several  other  cytokine  pathways  during  bacterial 
infection  and  it  is  assumed  that  similar  mechanisms  of  inflamma¬ 
tion  are  mediating  the  SIRS  to  HS.  Therefore,  using  a  conscious 
mouse  HS  model  we  measured  circulating  cytokine /soluble 
cytokine  receptor  (IL-1  a,  IL-1  /?,  IL-6,  IL-10,  TNF-a,  sIL-lRI, 
sIL-lRII,  sIL-6R,  sTNF-RI,  and  sTNF-RII),  liver  cytokine  and 
receptor  mRNA  accumulation  (IL-la,  IL- 1  /S,  IL-6,  IL-10,  TNF-a, 
IL-1RI,  IL-1RII,  IL-6R,  TNF-RI,  and  TNF-RII),  and  liver 
HSP70  mRNA  accumulation  in  wild-type  (WT)  and  TNFR  KO 
mice.  HSP70  mRNA  accumulation  levels  were  used  as  a  sensitive 
measure  of  stress  to  this  organ,  as  this  pathway  is  activated  by 
many  factors  that  are  inherent  in  our  HS  model  (e.g.,  heat  stress, 
dehydration,  oxidative  stress). 

The  mediators  involved  in  progression  of  the  HS  syndrome  and 
the  intricate  map  of  interactions  between  them  form  a  complex 
system  that  can  only  be  truly  understood  using  a  systems 
approach.  Although  several  computational  models  of  acute 
inflammation  exist  [19,20],  the  aim  of  this  study  was  to 
incorporate  a  level  of  mechanistic  detail  that  was  not  previously 
incorporated  into  these  models.  Therefore,  we  developed  a 
mathematical  model  that  integrates  relevant  biological  knowledge 
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with  our  novel  experimental  data  from  wild-type  mice  to  identify 
testable  hypotheses  that  will  delineate  the  molecular  mechanisms 
mediating  the  complex  etiology  of  the  heat-induced  SIRS.  This 
mechanistic  dynamic  model  describes  intra-  and  extracellular 
changes  (in  the  liver  and  in  the  plasma,  respectively)  in  cytokine 
signaling  pathways  under  HS  and  was  fitted  to  genomic  and 
proteomic  data  of  wild-type  mice  by  means  of  global  optimization 
techniques.  Model  validation  was  performed  using  a  completely 
different  set  of  data  from  TNFR  KO  mice  that  were  not  used  for 
calibration  purposes,  but  demonstrate  the  predictive  capabilities  of 
our  framework.  The  broader  applicability  of  the  developed  model 
in  the  context  of  acute  inflammation  was  assessed  by  comparing  its 
predictions  with  experimental  data  from  mice  treated  with  LPS 
[21].  The  purpose  of  this  study  was  to:  1)  assess  the  complex 
interaction  of  cytokines  in  the  liver  and  the  circulation  during  early 
progression  of  the  heat-induced  SIRS;  2)  gain  insight  into 
molecular  mechanism(s)  that  may  serve  as  future  therapeutic 
targets  for  HS  patients;  3)  analyze  the  correlation  of  organ  (liver) 
mRNA  accumulation  and  circulating  levels  of  cytokines;  4) 
determine  whether  the  TNFR  KO  responses  in  the  liver  and  the 
circulation  are  altered  due  to  differences  in  heating  and  cooling  or 
represent  a  direct  effect  of  the  absence  of  TNF  signaling. 
Moreover,  we  provide  a  unique  modeling  framework  that  supports 
identification  of  the  role  of  different  cytokine  signaling  pathways 
using  three  in  silico  experiments  that  can  be  used  to  guide  further 
in  vivo  experiments.  Specifically,  we  examined  the  response  to  the 
exposure  to  a  more  severe  heat  insult,  the  injection  of  a  dose  of 
LPS,  and  the  knockout  of  IL-1  OR. 

Results 

Mathematical  Model  of  the  Cytokines  Network 

The  primary  mechanisms  represented  in  this  model  are  the 
following:  activation  of  various  transcription  factors  (TFs)  by 
stimulation  via  a  set  of  external  and  internal  signals  triggered  by 
heat  stress;  regulation  of  cytokine  and  cytokine  receptor  gene 
transcription  involved  in  the  network  by  means  of  these  TFs; 
translation  of  mRNA  into  proteins;  transport  of  the  soluble 
proteins  to  the  pericellular  milieu;  binding  of  plasma  cytokines  to 
their  cognate  receptors;  and  signaling  back  to  the  TFs.  Details 
about  the  modeling  assumptions  are  given  in  the  Methods  section. 

The  elevated  Tc  in  response  to  heat  insult  is  one  of  several 
factors  that  is  thought  to  induce  organ  damage  and  increase  the 
concentration  of  denatured  proteins  (DP),  endotoxins  (lipopoly- 
saccharide,  LPS),  and  reactive  oxygen  species  (ROS)  that 
concomitantly  initiate  a  network  of  cytokine  responses.  Four  TFs 
were  assumed  to  be  the  primary  regulators  of  this  network, 
namely,  heat  shock  factor  1  (HSF-1),  nuclear  factor-K'B  (NF-kB), 
activator  protein  1  (AP-1),  and  signal  transducer  and  activator  of 
transcription  3  (STAT-3).  The  implicated  ligands  and  receptors 
were  classified  into  six  families  (heat  shock  protein  (HSP)70,  toll¬ 
like  receptor  (TLR)-4,  IL-1,  IL-6,  IL-10,  and  TNF)  and  their  inter- 
and  intracellular  activity  was  incorporated  into  the  model  as  six 
interconnected  modules  (Figure  1  and  Supporting  Information 
S2). 

The  model  is  described  by  a  system  of  65  ordinary  differential 
equations  (ODEs)  and  2 1 7  parameters  from  which  1 30  parameters 
were  fitted  to  experimental  data  corresponding  to  a  strain  of  WT 
mice  (B6129F2)  using  the  toolbox  SensSB  [22].  The  model  was 
built  based  on  the  hypothesis  that  the  liver  is  the  main  source  of 
circulating  cytokines.  In  order  to  test  this  hypothesis,  during  model 
fitting  we  gave  priority  to  the  liver  qPCR  data  (higher  weights  in 
the  cost  function);  therefore,  we  tried  to  achieve  the  best  possible 
fitting  for  die  liver  and  analyze  whether  it  was  possible  to 
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Figure  1 .  Cellular  network  of  interactions  amongst  HSP70,  TLR4,  IL-1,  IL-6,  IL-10,  and  TNF  families  induced  by  heat  stroke.  Prolonged 
exposure  to  a  high  ambient  temperature  increases  core  temperature  and  is  associated  with  organ  damage,  increase  of  denatured  proteins  (DP), 
reactive  oxygen  species  (ROS),  and  LPS.  As  a  consequence,  transcription  factors  (TF)  HSF-1,  NF-rcB,  AP-1,  and  STAT-3  are  activated  and  regulate  the 
transcription  of  a  set  of  cytokine  genes  in  the  cell  nucleus  represented  by  the  inner  grey  circle  (each  of  the  colored  boxes  contains  the  genes 
regulated  by  a  certain  TF:  HSF  in  yellow,  NF-jcB  in  green,  AP-1  in  blue,  and  STAT-3  in  dark  grey).  These  mRNAs  are  exported  to  the  cytoplasm  and 
translated  into  proteins  (not  represented  in  the  figure).  These  cytokines/soluble  cytokine  receptors  can  exit  the  cell  throughout  the  cell  membrane 
(outer  grey  circle)  and  be  released  into  the  blood  stream.  The  proteins  that  are  embedded  into  the  cell  membrane  represent  transmembrane 
receptors. 
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simultaneously  fit  the  circulating  cytokines  based  on  the  afore¬ 
mentioned  hypothesis.  For  the  mRNA  expression  of  the  WT 
animals,  the  fitting  lies  within  the  confidence  interval  of  almost 
every  experimental  point.  Only  the  values  of  mIL6  and  mILlO  for 
Tcjnax  are  considerably  smaller  than  those  predicted  by  the  model; 
disregarding  these  two  points,  the  average  of  the  absolute  values  of 
the  standardized  residuals  (e*  =  (_y,  —  JF/)/ cr)  is  smaller  than  one, 
meaning  that,  on  average,  the  error  of  the  fit  is  smaller  than  one 
standard  deviation  from  the  mean  [23] .  The  model  was  validated 
using  a  set  of  data  corresponding  to  TNFR  KO  mice  that  was  not 
used  for  model  calibration.  Note  that,  for  validation  purposes,  all 


parameters  were  fixed  to  the  best  values  obtained  in  the  calibration 
step  except  for  those  related  to  TNF-RI  and  II  transcription, 
which  were  set  to  zero.  A  file  with  the  value  of  the  parameters, 
reaction  rates,  and  model  equations  is  included  as  Supplementary 
Information  S3.  SensSB  files  needed  to  reproduce  the  results  are 
provided  as  Supporting  Information  S4.  Matlab  figures  containing 
the  estimated  time  course  for  all  the  states  are  available  as 
Supporting  Information  S5. 
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Experimental  Data,  Model  Fitting  and  Validation 

Liver  mRNA  accumulation  and  circulating  cytokines  and 
receptors  concentration  were  measured  at  four  sampling  points: 
1)  baseline  (prior  to  heat  exposure);  2)  maximum  core  temperature 
(7c,»mx  =  42.4  C),  also  referred  to  as  HS  collapse,  at  which  time 
the  mice  were  removed  from  the  heat;  3)  ~30  minutes  of  HS 
recovery  when  the  core  temperature  returns  to  baseline  (RTB;  first 
Tc  value  <36.0  C  during  cooling);  4)  ~3  hours  of  recovery  when 
mice  exhibit  hypothermia  depth  (lowest  Tc  value  with  cooling  rate 
of  0.01  C  *  min  during  recovery).  Liver  mRNA  accumulation  of 
AP-1  and  NF-kB  related  genes  is  shown  in  Table  1.  Liver  mRNA 
accumulation  of  HSP,  cytokines,  and  cytokine  receptors  is 
summarized  in  Table  2.  The  experimental  protocol  is  detailed  in 
the  Alethods  section. 

Mouse  core  temperature  profiles.  The  two  strains  of  mice 
used  in  this  study,  WT  and  TNFR  KO,  showed  significantly 
different  Tc  profiles  under  the  same  heat  stress  protocol  (Figure  2). 
TNFR  KO  mice  maintained  ~0.6°C  lower  Tc  than  WT  mice 
during  hyperthermia  but  experienced  similar  thermal  load  (a 
measure  of  heat  strain),  although  the  time  to  reach  Tc<max  was  ~30 
minutes  longer  than  the  WT  mice  (271  +  8%  vs.  240  +  8%, 

Table  1.  Liver  fold-change  in  mRNA  accumulation  during 
heat  stroke  recovery  in  WT  and  TNFR  KO  mice  for  NF-;cB  and 
AP-1  related  genes. 


Gene  Strain  T  c,max  RTB  Hypothermia 

NF-k-BI  B6129F2  0.750  0.872  1 . 1 4+a 

TNFR  KO  0.846  1.04  1 .29+Sa 

NF-kB2  B6129F2  1.37a  1.65  1.86 

TNFR  KO  1.55a  2.08  2.061 

RELA  B6129F2  1.08  1.08  1.461 

TNFR  KO  1.18  1.37  1.38 

RELB  B6129F2  1.55  2.02  2.03 

TNFR  KO  2.14s  2.14  3.16f 

C-REL  B6129F2  4.26  3.75  4.27 

TNFR  KO  4.54  3.69  4.32 

IkB  B6129F2  1.30  1.36  1.67 

TNFR  KO  1.05  1.42  1.52 

JUN  B6129F2  53.0a  109.9*  41.7 

TNFR  KO  40.8a  86.9*  93.01 

FOS  B6129F2  46.7  52.8  24.2a 

TNFR  KO  53.6  53.6  115.9ta 

Data  represent  fold-change  in  liver  mRNA  accumulation  relative  to  controls  at 
the  same  time  point.  Bold  fonts  represent  significantly  higher  than  time- 
matched  controls  (one-way  ANOVA,  P  <  0.05); 

^represents  significant  difference  from  T c,max) 

Represents  significant  difference  from  return-to-baseline; 

Represents  significant  difference  between  genotypes  within  each  gene.  TcmflX, 
maximum  core  temperature  (42.4°C);  TNFR  KO,  tumor  necrosis  factor  receptor 
knockout;  WT,  wild-type  strain  (B6129F2). 
doi:1 0.1 371 /journal. pone.0073393.t001 


Table  2.  Liver  fold-change  in  mRNA  accumulation  during 
heat  stroke  recovery  in  WT  and  TNFR  KO  mice  for  HSP, 
cytokines,  and  cytokine  receptors. 


Gene  Strain  T  Cmax  RTB  Hypothermia 

HSP70  B6129F2  657.9  780.41  1204.5ts 

TNFR  KO  565.8  783.21  1071.1ts 

alL-1  B6129F2  2.0a  3.0a  2.1 

TNFR  KO  1.4Sa  1.6ta  3.7+ 

[i  IL-1  B6129F2  2.6a  10.7ta  4.7 

TNFR  KO  2.2a  4.7ta  7.0+ 

IL-1RI  B6129F2  3.1  3.0  9.0ts 

TNFR  KO  3.0  2.9  4.7+ 

IL-1RII  B6129F2  0.6a  1.9a  9.4t§ 

TNFR  KO  0.5a  0.8ta  6.3ts 

IL-6  B6129F2  1.1 Sa  14.9ta  2.7Sa 

TNFR  KO  0.9Sa  6.7ta  6.4ta 

IL-6R  B6129F2  1.2  1.1  1.7s 

TNFR  KO  1.2  1.2  1.4+ 

gp130  B6129F2  1.4  1.0  1.5 

TNFR  KO  1.2  1.4  1.2 

IL-1 0  B6129F2  1.0  11.3ta  14.0* 

TNFR  KO  1.1  4.6+a  19.3+s 

aTNF-  B6129F2  0.4  2.0+  1.1  + 

TNFR  KO  0.5  0.7  0.8 

Data  represent  fold-change  in  liver  mRNA  accumulation  relative  to  controls  at 
the  same  time  point.  Bold  fonts  represent  significantly  higher  than  time- 
matched  controls  (one-way  ANOVA,  P  < 0.05); 

+ represents  significant  difference  from  T,m„„ ; 

Represents  significant  difference  from  return-to-baseline; 

Represents  significant  difference  between  genotypes  within  each  gene.  Tt.maxl 
maximum  core  temperature  (42.4  C);  TNFR  KO,  tumor  necrosis  factor  receptor 
knockout;  WT,  wild-type  strain  (B6129F2). 
doi:1 0.1 371/journal.pone.0073393.t002 

respectively;  ANOVA,  P<  0.05)  and  associated  with  a  significantly 
higher  level  of  dehydration  (11.8  +  0.3%  vs.  10.0  +  0.3%,  respec¬ 
tively;  ANOVA,  P<  0.05)  [24],  TNFR  KO  mice  also  showed  a 
significantly  faster  cooling  rate  than  WT  mice  from  Tc  max  to 
hypothermia  depth.  Despite  genotype  differences  in  heating  and 
cooling  profiles,  hypothermia  depth  was  virtually  identical 
between  genotypes  and  was  observed  ~3  hours  after  removal 
from  the  heat  stress  environment.  Taken  together,  these  data 
indicate  a  direct  effect  of  TNF  signaling  on  Tc  regulation  during 
heat  exposure  and  recovery  [24] . 

Liver  transcription  factor  mRNA  accumulation.  We 
measured  liver  mRNA  accumulation  of  all  the  NF-kB  family 
members,  namely,  NF-kB  l,  NF-kB2,  RELA,  RELB,  C-REL,  and 
IkB  (see  Table  1).  Liver  NF-kB  1,  NF-kB2,  RELA,  and  IkB 
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Figure  2.  Averaged  core  temperature  responses  for  WT  and  TNFR  KO  mice  during  heat  exposure  and  recovery  at  a  constant 
ambient  temperature.  Experimental  data  were  collected  at  four  different  points  based  on  mouse  core  temperature  responses  (baseline 
J'(.<36.0  C,  Tcmax  =  42.4'C,  return  to  baseline,  and  hypothermia  depth);  therefore,  the  average  was  computed  along  the  temperature  axis. 
doi:10.1371/journal.pone.0073393.g002 


mRNA  accumulation  was  not  significantly  or  little  altered  (less 
than  two-fold)  at  any  time  point  in  either  genotype.  C-REL  was 
upregulated  ~4-fold  at  all  sampling  time  points  after  the  onset  of 
HS  (Tc,max,  RTB,  and  hypothermia  depth)  and  RELB  was  ~2- 
fold  increased  at  RTB  and  hypothermia  with  less  mRNA 
accumulation  observed  at  Tcmax  for  the  WT  mice.  Although  we 
did  not  direcdy  measure  NF-rcB  nuclear  localization,  some  studies 
suggest  that  increased  RELB  transcription  and  translation  is  a 
direct  measure  of  this  activity  [25] .  Overall,  we  did  not  detect  any 
genotype  differences  in  liver  mRNA  accumulation  of  the  NF-rcB 
family  members. 

The  most  robust  changes  between  genotypes  in  the  liver  were 
observed  in  JUN  and  FOS  mRNA  accumulation  profiles  (Table  1). 
While  both  TFs  were  significantly  upregulated  (fold-change  ~55.3 
and  ~38.5,  respectively)  at  TCymax,  JUN  mRNA  accumulation  was 
significantly  lower  at  Tc<max  in  TNFR  KO  compared  to  WT  mice. 
Liver  JUN  mRNA  accumulation  reached  peak  levels  at  RTB  in 
the  WT  mice  (fold-change  — 123)  whereas  TNFR  KO  mice 
showed  delayed  mRNA  accumulation  with  maximum  values 
observed  at  hypothermia  (fold-change  ~92).  Interestingly,  JUN 
and  FOS  mRNA  accumulation  at  baseline  was  slightly  but 
significantly  lower  (fold-change  ~0.8)  in  the  liver  of  TNFR  KO 
compared  to  WT  mice  suggesting  homeostatic  regulation  of  this 
response  by  TNF. 

HSPs.  Baseline  liver  FISP70  mRNA  accumulation  did  not 
differ  between  genotypes  and  values  observed  in  non-heated 
controls  were  similar  across  all  time  points  (data  not  shown).  At 
Tcjnax,  both  genotypes  showed  ~600-fold  increase  in  liver  HSP70 
mRNA  accumulation  that  was  significantly  higher  than  controls, 
similar  between  genotypes,  and  sustained  through  RTB  in  both 
groups  (Table  2;  one-way  ANOVA,  P<  0.001).  Peak  mRNA 
accumulation  of  HSP70  was  observed  at  hypothermia  with  no 
significant  difference  between  genotypes  (WT:  1364-fold; 
TNFRKO:  1034-fold;  Table  2). 

Figure  3  shows  model  prediction  based  on  WT  data  versus 
experimental  data  from  TNFR  KO  mice  for  liver  mRNA 
accumulation  of  HSP70.  In  accordance  with  literature  [26]  as 
well  as  our  experimental  data,  the  model  accurately  predicts 
increased  liver  HSP70  mRNA  accumulation  shortly  after  HS 
collapse  (~4— 4.5  hours),  with  peak  mRNA  accumulation  (>1000 
fold)  occurring  at  ~7  hours,  when  mice  are  hypothermic.  As 
illustrated  in  Figure  3,  our  model  accurately  fits  the  WT  data  used 


for  calibration  (solid  purple  line)  and,  more  importantly,  predicts  a 
slightly  delayed  and  attenuated  response  of  the  TNFR  KO  strain 
(dashed  green  line).  Although  our  experimental  data  from  TNFR 
KO  mice  suggested  attenuated  liver  HSP70  mRNA  accumulation, 
this  trend  was  not  statistically  significant.  However,  one  might 
expect  such  a  response  to  occur  in  the  TNFR  KO  mice  since  the 
higher  level  of  dehydration  experienced  by  these  mice  would 
generate  greater  ROS  production  and  consequent  inhibition  of 
HSF-1. 

IL-1  family.  Figure  4  represents  model  predictions  (solid  and 
dashed  lines)  and  experimental  data  (markers)  for  IL- 1  a.  and  IL- 1  /? 
at  the  four  sampling  points  (baseline,  T c.max,  RTB,  and 
hypothermia).  The  left  panels  (4A  and  4C)  correspond  to  model 
fitting  using  WT  data,  while  right  panels  (4B  and  4D)  illustrate 
model  validation  using  TNFR  KO  data  that  were  not  used  for 
calibration.  The  top  panels  (4A  and  4B)  represent  the  fold  change 
for  liver  mRNA  accumulation  (mIL- 1  a,  mlL- 1  /J)  and  the  bottom 
panels  (4C  and  4D)  correspond  to  the  concentration  of  soluble  (i i.e ., 
circulating)  cytokines  (sIL-la,  sIL-ljS). 

IL-1  a,  IL-1  /?,  and  IL-IRa  promoters  contain  NF-fcB  and  AP-1 
regulatory  elements  [27,28].  However,  IL-1  a  and  IL-1  /?  exhibited 
differential  mRNA  accumulation  in  the  liver  of  HS  mice,  with  IL- 
1  fi  fold  change  significantly  higher  than  IL- 1  a  at  all  time  points 
(see  Table  2  and  Figures  4A-B).  Notably,  liver  IL-1  a  and  IL-1 /I 
mRNA  accumulation  was  significantly  delayed  in  TNFR  KO 
mice;  these  experimental  data  were  captured  by  the  model 
(Figure  4B).  That  is,  our  model  predicted  an  earlier  peak  of  mlL- 
1  a  and  ft  mRNA  accumulation  in  the  WT  at  RTB  (third  sampling 
point)  at  ~5  hours  while  TNFR  KO  mice  presented  maximal 
mRNA  accumulation  between  RTB  and  hypothermia  (third  and 
fourth  sampling  point)  at  ~6  hours.  IL-1RI  mRNA  accumulation 
was  similarly  elevated  in  both  genotypes  throughout  recovery  (~3- 
fold  increase),  but  showed  peak  mRNA  accumulation  at  hypo¬ 
thermia  (~8-fold  for  WT  and  ~4-fold  for  TNFR  KO)  compared 
to  the  earlier  time  points.  IL-1RII  mRNA  accumulation  was 
similar  between  groups  with  a  significant  increase  in  heated 
animals  compared  to  controls  at  hypothermia  only  (fold  change 
~6). 

The  concentration  of  sIL- 1  /?  (i.e.,  circulating  protein)  was  ~  10- 
fold  higher  than  sIL- 1  a  at  all  time  points  regardless  of  genotypes 
(Figures  4C-D).  This  difference  in  plasma  levels  might  be  due  to 
the  different  signaling  mechanisms  of  IL- 1  a  and  IL- 1  /?  that  induce 
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Figure  3.  Model  simulation  (solid  and  dashed  lines)  versus  experimental  data  (markers)  for  the  mRNA  accumulation  of  HSP70 
under  heat  stroke.  The  model  accurately  fits  the  WT  data;  moreover,  it  predicts  a  slightly  delayed  and  attenuated  response  for  the  TNFR  KO  mice  in 
agreement  with  the  experimental  data  (data  used  for  validation  only). 
doi:10.1371/journal.pone.0073393.g003 

expression  of  their  own  genes  in  an  autocrine  and  paracrine 
manner,  respectively  [29] .  In  addition,  HS  had  no  significant  effect 
on  sIL-la  sIL-l/f,  sIL-lRI  or  sIL-lRII  levels  in  either  genotype. 

The  dissociation  between  liver  IL-1  family  mRNA  accumulation 
and  circulating  protein  levels,  reproduced  by  the  computational 
model,  suggests  that  the  proteins  generated  in  the  liver  are  not 
extensively  released  into  the  blood  stream  or  present  a  long  delay 
not  captured  within  the  time  frame  of  our  experimental  protocol. 

IL-6  family.  Figure  5  shows  the  dynamics  of  IL-6  and  IL-6R 
with  a  structure  analogous  to  Figure  4  (left  panels  for  WT,  right 
panels  for  TNFR  KO;  top  panels  for  mRNA  accumulation  fold 
change,  bottom  panels  for  circulating  concentration).  WT  liver 
mIL-6R  was  not  significantly  affected  by  HS  (dashed  line);  in 


contrast,  mIL-6  was  undetectable  at  baseline,  but  ~  13-fold  higher 
than  the  limit  of  detection  at  RTB,  which  represented  the  peak 
time  of  mRNA  accumulation  of  this  cytokine.  The  model 
accurately  fits  the  experimental  data  for  both  IL-6R  and  IL-6 
(Figure  5A)  and  model  validation  showed  that  the  model 
accurately  predicts  the  TNFR  KO  response.  That  is,  the  model 
predicts  a  delay  with  respect  to  WT  data  with  maximum  mRNA 
accumulation  between  RTB  and  hypothermia  (~6  hours)  in 
contrast  to  the  maximum  at  RTB  experienced  by  WT  animals  (~5 
hours,  Figures  5A-B).  The  model  predictions  for  sIL-6  are 
reasonable  with  peak  levels  expected  at  ~5-6  hours  (Figures  5C— 
D).  On  the  contrary,  sIL-6R  responses  were  dissociated  from  the 
time  course  and  magnitude  of  changes  observed  in  mIL-6R 
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Figure  4.  Model  simulation  (solid  and  dashed  lines)  versus  experimental  data  (markers)  for  the  IL-1  family.  (A)  represents  the  mRNA 
accumulation  of  IL-1  a  and  IL-1 /I  for  WT,  (B)  the  mRNA  accumulation  of  IL-1  a  and  IL-1  /I  for  TNFR  KO  (data  used  for  validation  only),  (C)  the  plasma 
concentration  of  IL-1  a  and  IL-1  fi  for  WT,  and,  (D)  the  plasma  concentration  of  IL-1  a  and  IL-1 /I  for  TNFR  KO  (data  used  for  validation  only). 
doi:1 0.1 371 /journal. pone.0073393.g004 
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Figure  5.  Model  simulation  (solid  and  dashed  lines)  versus  experimental  data  (markers)  for  the  IL-6  family.  (A)  represents  the  mRNA 
accumulation  of  IL-6  and  IL-6R  for  WT,  (B)  the  mRNA  accumulation  of  IL-6  and  IL-6R  for  TNFR  KO  (data  used  for  validation  only),  (C)  the  plasma 
concentration  of  IL-6  and  slL-6R  for  WT,  and  (D)  the  plasma  concentration  of  IL-6  and  slL-6R  for  TNFR  KO  (data  used  for  validation  only). 
doi:1 0.1 371  /journal. pone.0073393.g005 


indicating  that  the  liver  is  not  an  important  source  of  this 
circulating  protein.  As  a  consequence,  the  model  is  not  able  to 
accurately  predict  the  experimental  changes  in  sIL-6R  for  either 
WT  or  TNFR  KO. 

IL-10  family.  In  our  experimental  mouse  study,  liver  mIL- 1 0 
was  undetectable  at  baseline  and  peaked  at  hypothermia  showing 
~  20-fold  increase  above  the  limit  of  detection  observed  in  both 
WT  and  TNFR  KO  mice  (see  Figures  6A-B);  this  response  was 
accurately  captured  by  our  model.  We  observed  dissociation 
between  mIL-10  and  sIL-10  levels  that  the  mathematical  model 
was  not  able  to  capture  (see  Figures  6C-D).  That  is,  liver  mlL- 
10  mRNA  accumulation  differed  between  strains  at  RTB  with 
significantly  higher  levels  observed  in  WT  compared  to  TNFR 
KO  mice  (Table  2  and  Figures  6A-B).  On  the  other  hand,  we 
have  observed  significantly  lower  sIL- 1 0  levels  in  WT  compared  to 
TNFR  KO  mice  at  hypothermia  ( — 150  and  ~500  pg/ml, 
respectively;  [24]),  which  was  not  accounted  for  by  our  model 
(Figures  6C-D).  This  disconnect  between  predicted  and  experi¬ 
mental  sIL-10  levels  might  indicate  that  the  mRNA  accumulated 
in  the  liver  is  not  consistently  translated  into  proteins  and/or 
transported  to  the  plasma,  as  noted  previously  for  some  of  the  IL- 1 
and  IL-6  family  members. 

TNF-a  family.  Liver  mTNF-a  did  not  show  significant 
change  in  response  to  HS  in  either  genotype  (Table  2).  TNFR 
KO  mice  showed  ~6-fold  higher  sTNF-a  levels  than  WT  mice  at 
baseline  and  all  time  points  of  HS  recovery  in  our  model,  with  no 
effect  of  heat  exposure  on  this  cytokine  (Figure  7,  [24]).  Our  model 
slightly  overestimates  the  sTNF-a  concentration  observed  for  the 
WT  mice  and  predicted  a  sTNF-a  profile  for  the  TNFR  KO  mice 
significandy  higher  than  the  one  observed  experimentally 
(Figure  7A).  Liver  mTNF-RI  and  mTNF-RII  were  not  signifi¬ 
cantly  altered  by  HS;  yet,  sTNF-RI  and  sTNF-RII  levels  showed  a 
significant  increase  during  HS  recovery  (up  to  ~5  fold-change  for 
sTNF-RI  at  TCimax;  see  Figure  7B).  The  failure  of  our  model  to 
capture  this  robust  change  in  sTNFRI  and  II  indicates  that  the 
liver  profile  was  not  accurately  reflecting  circulating  changes  in 


protein;  as  such,  the  liver  did  not  appear  to  be  the  main  source  of 
these  soluble  receptors.  As  expected,  plasma  levels  of  the  sTNF-RI 
and  RII  proteins  were  undetectable  at  all  time  points  in  TNFR 
KO  mice. 

Sensitivity  and  Correlation  Analysis 

Despite  the  use  of  global  optimization  techniques,  we  found 
several  sets  of  parameters  that  accurately  fit  the  data,  with 
differences  oil  the  order  of  the  precision  of  the  experimental  data. 
Some  of  the  sets  lead  to  indistinguishable  model  predictions  that 
can  be  explained  by  the  correlation  between  parameters,  the  low 
sensitivity  of  some  of  them,  and/ or  different  alternative  pathways 
that  the  available  data  are  not  able  to  discriminate.  The  set  with 
the  smallest  objective  function  value  was  selected  as  the  optimal  set 
and  used  for  the  simulations  depicted  in  the  figures  and  further 
analyses. 

Local  relative  sensitivities  revealed  that,  for  the  optimal 
parameter  set,  30%  of  the  parameters  account  for  less  than  2% 
of  the  information  while  the  12%  most  influential  account  for 
more  than  50%  of  the  total  sensitivity.  The  most  influential 
parameters  are  those  related  to  HSP  transcription,  LPS  action 
(initial  LPS  content  in  the  guts  and  signaling  through  TLR4),  and 
IL-6  and  IL- 1  /J  transcriptional  activation  mediated  by  AP- 1  and 
NF-kB. 

Figure  8  depicts  the  relative  sensitivity  indices  (SI)  and  the 
correlation  matrix  for  the  parameters  of  the  Hill  functions  involved 
in  AP-1  target  genes  transcriptional  activation.  SI  show  that  IL- 
lRa-related  parameters  ( kjnIL\rci_AP\ ,  kkjnIL\ra_APl,  and 
nj?iIL\ra_APl)  are  the  least  influential  whereas  IL-l/f-related 
parameters  ( k_mIL\b_AP\ ,  kkjnIL\b_AP\,  and  njnIL\b_APV) 
are  the  most  important;  this  is  in  agreement  with  the  experimental 
data,  which  indicated  that  IL- 1  fl  was  the  main  member  of  the  IL- 1 
family  with  actions  during  HS  recovery.  The  correlation  matrix 
shows  strong  correlation  among  the  three  parameters  involved  in 
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Figure  6.  Model  simulation  (solid  and  dashed  lines)  versus  experimental  data  (markers)  for  IL-1 0.  (A)  represents  the  mRNA  accumulation 
of  IL-1 0  for  WT,  (B)  the  mRNA  accumulation  of  IL-1 0  for  TNFR  KO  (data  used  for  validation  only),  (C)  the  plasma  concentration  of  IL-1 0  for  WT,  and  (D) 
the  plasma  concentration  of  IL-1 0  for  TNFR  KO  (data  used  for  validation  only). 
doi:1 0.1 371  /journal. pone.0073393.g006 


each  of  the  Hill  equations  (see  Eq.  6),  which  may  be  the  source  of 
some  of  the  identifiability  issues. 

Strong  correlations  were  also  found  among  the  parameters  that 
regulate  AP- 1  activation  by  ROS  and  LPS.  This  is  due  to  the  lack 
of  ROS  and  LPS  experimental  data,  which  makes  the  distinction 
between  the  two  main  mechanisms  for  AP-1  activation  very 
difficult. 

In  silico  Experiments 

To  test  the  strength  and  utility  of  our  model  for  analysis  of  the 
heat-induced  SIRS,  we  conducted  a  series  of  in  silico  experiments 
testing  the  effect  of  higher  heat  stress  temperature,  and  knockout 
of  genes  that  have  been  implicated  in  the  inflammatory  response 
(i.e.,  IL-1  OR).  Moreover,  an  LPS  injection  experiment  similar  to 
the  one  performed  in  [2 1]  was  simulated  and  the  predicted  values 


for  liver  mRNA  accumulation  and  plasma  concentration  for 
several  cytokines  were  qualitatively  compared  with  the  mice  data. 

Severity  of  the  Heat  Insult 

Several  studies  show  divergent  results  with  respect  to  the 
circulating  TNF-a  response  to  HS.  Hammami  et  al.  [9]  were 
unable  to  detect  elevated  circulating  TNF-a  concentrations  in 
2 1  HS  patients  at  the  time  of  clinical  presentation  or  shortly  after 
cooling.  On  the  contrary,  Bouchama  et  al.  [30]  observed  high 
circulating  TNF-a  levels  at  clinical  admission  in  17  HS  patients.  It 
is  unclear  if  incongruity  between  studies  is  a  complication  of 
differences  in  experimental  procedures,  extent  of  heat  injury  or 
other  unidentified  factors.  Several  experiments  have  shown  that 
HSP70  (induced  by  heat  or  oral  glutamine)  inhibits  TNF-a  [31- 
34] .  Moreover,  Shell  et  al.  [35]  have  observed  that  heat  shock 
inhibits  NF-kB  activation  in  a  dose-  and  time-dependent  manner. 
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Figure  7.  Model  simulation  (solid  and  dashed  lines)  versus  experimental  data  (markers)  for  the  TNF  family.  (A)  represents  the  plasma 
concentration  of  TNF-a  for  WT  and  TNFR  KO,  and  (B)  the  plasma  concentration  of  the  soluble  isoforms  of  receptors  TNF-RI  and  TNF-RII  for  WT. 
doi:10.1371/journal.pone.0073393.g007 
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Figure  8.  Sensitivity  indices  and  correlation  matrix  for  the  Hill  function  parameters  of  the  AP-1  target  genes.  SI  show  that  IL-IRa 
related  parameters  are  the  least  influential  whereas  IL-1  /f  related  are  the  most  important.  The  correlation  matrix  shows  strong  correlations  among  the 
three  parameters  involved  in  each  of  the  Hill  equations. 
doi:1 0.1 371  /journal. pone.0073393.g008 


An  in  silico  experiment  using  a  more  severe  heat  insult  as  an  input 
to  the  model  resulted  in  increased  TNF-a  mRNA  accumulation 
levels  (Figure  9).  In  response  to  Tc  max  of  42.4°C,  which  simulates 
the  heat  severity  experienced  by  WT  and  TNFR  KO  mice  in  our 
in  vivo  study,  we  observed  decreased  TNFa  mRNA  accumulation 
throughout  HS  recovery  (Figure  9).  Interestingly,  no  significant 
changes  in  TNFa  mRNA  accumulation  were  observed  during  the 
initial  8  hours  of  experimentation,  which  fits  our  experimental 
data.  Conversely,  increasing  the  TCtmax  to  45.0°C  induced  ~4  fold 
increase  in  TNFa  mRNA  accumulation  within  ~16  hours  of 
recovery.  This  in  silico  experiment  supports  the  hypothesis  in  [35] 
that  high  or  low  concentrations  of  TNF-a  are  a  function  of  the 
duration  and  severity  of  the  heat  insult  that  might  also  influence 
the  severity  of  the  damage  to  the  gut  epithelial  barrier  and 
therefore  the  release  of  LPS. 


IL10-R  KO.  The  anti-inflammatory  actions  of  IL-10  are 
thought  to  be  mediated  through  the  suppressor  of  cytokine 
signaling  3  (SOCS-3)  pathway.  An  in  silico  KO  of  the  IL-1  OR 
showed  that  SOCS-3  signaling  was  significantly  decreased 
compared  to  the  response  observed  in  WT  mice  (Figure  10).  That 
is,  WT  mice  showed  peak  accumulation  of  SOCS-3  mRNA  at  ~9 
hours;  although  the  time  course  of  this  response  was  similar  in  the 
IL-1  OR  KO,  SOCS-3  signaling  was  significantly  attenuated 
compared  to  the  WT  condition  (Figure  10).  These  data  suggest 
that  SOCS-3  signaling  is  a  downstream  target  of  IL-10  signaling 
that  may  regulate  the  inflammatory  response  during  recovery. 

LPS  injection.  In  order  to  assess  the  validity  of  our  model  in  a 
broader  context  of  systemic  inflammation  we  have  simulated  the 
experiment  described  in  [21]  where  WT  (C57BL/10)  mice  were 
treated  with  bacterial  LPS  and  liver  mRNA  levels  for  TNF-a  and 
ILl/J,  among  other  cytokines,  were  measured.  The  model 


Figure  9.  Predicted  TNF-a  mRNA  accumulation  for  WT  mice  under  different  temperature  profiles.  An  in  silico  experiment  using  a  more 

severe  heat  insult  as  an  input  to  the  model  resulted  in  increased  TNF-a  mRNA  accumulation  levels. 

doi:10.1371/journal.pone.0073393.g009 
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Figure  10.  Predicted  SOCS-3  mRNA  accumulation  for  WT  and  IL-10R  KO.  An  in  silico  KO  of  the  IL-10R  shows  decreased  levels  of  SOCS-3 
which  may  result  in  an  enhanced  pro-inflammatory  response. 
doi:10.1371/journal.pone.0073393.g010 

reproduces  the  experimental  results  showing  a  fast  inflammatory 
response  with  rapid  increase  of  TNF-a  and  IL 1  /?  liver  mRNA 
(Figure  11).  The  level  of  TNF-a  experiences  a  significantly  higher 
increase  than  under  heat  stroke  experiment  suggesting  that  LPS 
might  not  be  playing  a  crucial  role  in  the  SIRS  following  the  heat 
stroke  protocol  considered  in  this  work.  However,  strain  differ¬ 
ences  between  studies  and  the  use  of  only  one  dose  of  LPS 
complicate  the  interpretation.  Future  studies  will  need  to  be 
designed  to  determine  the  role  of  LPS  as  an  initiating  factor  of  the 
SIRS  in  our  heat  stroke  model. 

Discussion 

This  is  the  first  study  to  devise  a  complex  mathematical  model 
of  the  intra-  and  extracellular  cytokine  signaling  pathway 
interactions  that  are  stimulated  in  response  to  HS.  Literature 
was  thoroughly  reviewed  and  the  postulated  model  integrates 
relevant  biological  knowledge  and  novel  experimental  data.  Tissue 
damage,  DP,  LPS,  and  ROS,  stimulated  upon  Tc  increase,  are 
assumed  to  initiate  a  network  of  HSP  and  cytokine  responses 


involving  HSP70  and  four  families  of  cytokines,  namely  IL- 1,  IL-6, 
IL-10,  and  TNF.  Four  transcription  factors  (HSF-1,  NF-kB,  AP-1, 
STAT-3)  were  considered  to  be  the  primary  regulators  of  this 
network.  The  biological  assumptions  made  for  building  this  model 
are  supported  by  solid  scientific  evidence.  However,  it  has  to  be 
noted  that  the  relation  among  all  the  species  is  not  exacdy  known 
and  further  experimentation  might  invalidate  some  of  the 
assumptions.  Moreover,  model  fitting  and  sensitivity  analysis 
provide  important  information  about  the  strength  of  the  postulated 
interactions  giving  flexibility  to  the  model  for  supporting  or 
questioning  them. 

Hill  functions  and  mass  action  kinetics  were  used  to  build  an 
ODEs  model  that  was  calibrated  using  a  set  of  genomic  (liver 
mRNA  accumulation)  and  proteomic  (circulating  cytokines  and 
receptors)  data  by  means  of  global  optimization.  Despite  different 
factors,  as  blood  volume  changes  due  to  dehydration,  could 
influence  the  value  of  the  model  parameters,  these  effects  were 
considered  negligible  and  all  the  parameters  were  considered 
invariant  during  the  experiment.  The  model  was  validated  using 
an  independent  set  of  data  (not  used  for  calibration)  from  TNFR 


Figure  11.  Predicted  TNF-a  and  IL1  />  mRNA  accumulation  after  an  LPS  injection.  A  simulation  of  an  LPS  injection  in  WT  mice  shows  an 
immediate  inflammatory  response  with  rapid  increase  of  TNF-a  and  IL1/1  liver  mRNA. 
doi:10.1371/journal.pone.0073393.g01 1 
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KO  mice.  For  both  WT  and  TNFR  KO,  the  model  fits  the  liver 
qPCR  data  within  the  precision  of  the  experimental  data; 
however,  for  some  of  the  cytokine  receptors  (e.g.  sILGR,  sTNF- 
RI,  sTNF-RII)  there  is  no  combination  of  parameters  that 
accurately  fit  the  liver  and  plasma  data  simultaneously  since  the 
circulating  levels  are  elevated  while  the  liver  mRNA  remains 
unchanged.  The  deficiencies  in  the  model  fit  with  respect  to  these 
soluble  cytokine  receptors  led  to  this  relevant  conclusion: 
molecules  newly  transcribed  and  translated  in  the  liver  are 
unlikely  to  be  the  primary  sources  of  sIL6R,  sTNF-RI,  and 
sTNF-RII  in  the  circulation  during  HS.  Since  the  model  was  built 
based  on  this  hypothesis,  diis  fitting  mismatch  can  only  be 
overcome  by  changing  the  model  structure  and  including  other 
sources  of  circulating  cytokines.  Although  changes  of  clearance 
and/ or  degradation  rates  during  HS  could  eventually  explain  this 
divergence,  the  most  plausible  hypothesis  is  that  other  organs/cell 
types  (circulating  macrophages,  splenocytes),  which  we  plan  to 
integrate  into  a  more  comprehensive  model  in  the  future,  are  also 
contributing  to  the  circulating  cytokine  milieu.  As  we  move 
forward  with  this  modeling  effort,  these  insights  will  be  critically 
important  for  us  to  identify  the  sources  of  circulating  cytokines  and 
potential  therapeutic  targets  that  can  be  tested  physiologically  in 
our  mouse  model.  Moreover,  these  findings  motivate  the 
development  of  a  multi-modular  model  by  systematic  incorpora¬ 
tion  of  other  relevant  peripheral  organs,  such  as  the  spleen  and  the 
brain,  and  ultimately  the  entire  body. 

For  some  of  the  species  (e.g.  IL-1  /?,  IL-1  a,  IL-10),  changes  in 
liver  mRNA  accumulation  are  not  reflected  in  the  circulating 
protein  levels.  Although  our  computational  model  is  able  to 
capture  these  dynamics,  the  low  values  estimated  for  the 
translation/ transport  parameters  associated  with  these  cytokines 
suggest  that  the  mRNA  is  not  consistently  translated  into  proteins, 
the  proteins  generated  in  the  liver  are  not  extensively  released  into 
the  blood  stream,  and/or  they  present  a  long  delay  not  captured 
within  the  time  frame  of  our  experimental  protocol.  To 
discriminate  between  these  scenarios,  studies  on  protein  transla¬ 
tion  including  liver  protein  concentration  assays  and  protein 
tagging  should  be  performed. 

NF-kB  remains  practically  unchanged  in  the  liver  of  WT  and 
TNFR  KO  mice  under  the  experimental  HS  protocol  used  in  this 
study.  These  findings  are  in  agreement  with  those  of  [36]  that 
showed  that  NF-kB  DNA  binding  activity  did  not  change  in 
response  to  heat  stress  in  rats;  yet,  they  are  somewhat  surprising 
since  NF-kB  has  been  implicated  as  the  main  TF  that  regulates 
cytokine  mRNA  accumulation  during  inflammation  [37].  More¬ 
over,  both  TNF-a  liver  mRNA  accumulation  and  plasma 
concentrations  (activated  by  NF-kB)  were  not  affected  by  the 
heat  insult  in  contrast  to  other  studies  that  reported  elevated 
concentrations  of  this  cytokine  after  HS  [30].  Although  these 
discrepancies  between  studies  may  be  partly  explained  by  species 
differences  in  HS  responses,  a  more  plausible  explanation  is  that 
the  highly  elevated  mRNA  accumulation  of  HSP70  inhibited  the 
translocation  of  NF-kB  to  the  nucleus  in  a  dose-  and  time- 
dependent  manner  [35].  We  hypothesize  that  activation  of  NF-kB 
and  increased  concentration  of  TNF-a  depend  on  the  duration 
and  severity  of  the  heat  insult.  An  in  silico  experiment  with  a  more 
severe  heat  insult  illustrates  that  our  model  supports  this 
hypothesis  with  a  molecular  mechanism  (inhibition  of  NF-kB  by 
HSP70)  capable  of  explaining  the  different  behaviors,  encompass¬ 
ing  the  apparently  contradictory  results  from  the  literature. 

AP- 1  is  revealed  here  as  the  main  transcription  factor  activated 
in  the  liver  during  HS  recovery  and  it  appears  to  be  responsible  for 
the  high  expression  of  many  pro-inflammatory  cytokines  (i.e.  IL- 
1/1,  IL-6,  and  IL-10)  in  spite  of  the  low  activation  of  NF-kB.  The 


fact  that  JUN  and  FOS  liver  mRNA  accumulation  is  lower  in 
TNFR  KO  mice  at  baseline  indicates  that  the  differences  in 
TNFR  KO  responses  are  not  only  due  to  differences  in  heating 
and  cooling  but  a  direct  effect  of  the  lack  of  TNF  signaling.  It  is 
intriguing  to  speculate  that  lower  JUN  and  FOS  mRNA 
accumulation  may  be  a  direct  consequence  of  the  absence  of  a 
TNF  autocrine  cascade  effect  on  AP- 1  activation,  which  has  been 
shown  to  be  triggered  in  hepatocytes  in  response  to  TNF  [38] .  IL- 
1  /i,  IL-1RI,  IL-6,  and  IL-6R  also  presented  lower  mRNA 
accumulation  in  normothermic  conditions  for  TNFR  KO,  which 
can  be  explained  by  the  lower  activation  of  AP- 1 .  However,  we 
reject  the  hypothesis  that  the  changes  in  liver  gene  expression  in 
our  model  were  only  related  to  the  TNF  autocrine  cascade  since 
we  failed  to  detect  changes  in  liver  TNF  mRNA  accumulation  in 
WT  mice  at  any  time  point  of  HS  recovery.  Rather,  our  data 
suggest  that  delayed  liver  gene  expression  was  due  to  a 
combination  of  factors  that  include:  1)  lower  initial  (baseline) 
expression;  2)  indirect  effect  of  other  genes  that  activate  liver  AP- 1 
that  also  showed  a  delayed  response  {i.e.  IL-6,  IL-1);  3)  lower  core 
temperature  of  TNFR  KO  mice  during  heat  exposure.  Unfortu¬ 
nately,  it  was  not  possible  to  normalize  the  heating  responses 
between  genotypes  since  the  absence  of  TNF  signaling  induced  a 
downward  shift  in  the  temperature  set  point,  thus  causing  the  KO 
mice  to  regulate  Tc  and  metabolic  rate  differentially  from  WT 
mice  [24], 

AP-1  may  deactivate  the  heat  shock  response  during  stress 
recovery  by  hyperphosphorylating  HSF-1  to  inhibit  its  function 
[39,40];  hence,  robust  activation  of  AP-1  (e.g.  due  to  an  excess  of 
LPS)  could  lead  to  considerably  decreased  HSP70  gene  expression 
in  response  to  HS.  Since  HSP70  is  concomitantly  inhibiting  the 
activation  of  NF-kB,  HSP70  downregulation  could  exacerbate  the 
situation  leading  to  a  lethal  concentration  of  proinflammatory  NF- 
;cB  dependent  genes,  such  as  TNF-a.  Therefore,  AP-1  and  its 
related  cofactor  genes  (Jun  and  Fos  families)  are  revealed  here  as 
promising  targets  for  future  therapeutic  intervention  to  accelerate 
HS  recovery.  In  contrast,  HSPs  exert  an  anti-inflammatory  effect 
by  inhibiting  the  translocation  of  NF-kB  to  the  nucleus;  therefore, 
it  is  expected  that  drug-induced  increases  in  HSP70  gene 
expression  prior  to  or  following  HS  collapse  may  attenuate  the 
heat-induced  SIRS  in  the  liver  and  perhaps  other  organ  systems. 

An  important  outcome  of  our  predictive  model  is  the 
identification  of  a  differential  time  course  of  circulating  cytokine 
responses  in  TNFR  KO  mice  compared  to  their  WT  controls. 
This  is  an  important  aspect  of  our  model  as  it  identified  unique 
time  windows  that  should  be  considered  in  future  analyses  of  the 
circulating  biomarkers  that  may  be  mediating  damage.  Further¬ 
more,  it  suggests  that  alterations  of  the  cytokine  balance  (i.e.,  WT 
vs.  TNF  KO)  may  skew  the  cytokine  milieu  towards  a  pro-  or  anti¬ 
inflammatory  phenotype  that  alters  the  time  course  of  progression 
of  the  SIRS.  Future  studies  in  our  mouse  HS  model  that  target 
different  time  points  or  core  temperatures  for  cytokine  analysis  will 
be  instrumental  in  defining  the  rapid  changes  in  cytokine 
production  that  mediate  changes  in  our  model  and  may  have 
been  missed  in  our  current  analysis. 

The  uneven  distribution  of  the  parameter  sensitivities  found  in 
this  study  reinforces  Gutenkunst  et  al.  [41]  who  concluded,  after 
testing  several  systems  biology  models,  that  “sloppy”  spectra  of 
parameter  sensitivities,  i.e.  with  eigenvalues  roughly  evenly 
distributed  over  many  decades,  are  universal  in  systems  biology 
models.  This  property  may  explain  the  difficulty  of  extracting 
precise  parameter  estimates  from  collective  fits  and  reinforces  the 
need  for  establishing  a  parameter  ranking.  For  that  reason,  in  this 
work  the  parameter  estimation  was  done  in  two  stages;  first  we 
focused  on  the  most  influential  group  of  parameters,  whereas  the 
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less  important  group  was  fitted  in  a  second  stage.  Despite  the 
sequential  parameter  identification  and  the  use  of  global 
optimization  techniques,  the  identifiability  analysis  revealed  a 
number  of  difficulties  with  estimating  a  unique  value  for  the 
parameters.  Although  model  reduction  could  be  attempted  at  this 
stage,  we  preferred  to  analyze  the  detailed  mechanistic  model  and 
exploit  the  identifiability  deficiencies  encountered  for  planning 
future  experiments  aiming  to  obtain  a  complete  picture  of  the 
SIRS  ensuing  during  HS  recovery.  In  particular,  further 
experimentation  that  assesses  changes  in  circulating  LPS  concen¬ 
trations  or  the  effect  of  neutralizing  its  effects  (with  antibiotics  or 
gut  sterilization  treatments)  will  be  important  to  improve 
identifiability  and  discriminate  among  ROS  and  LPS  overlapping 
mechanisms.  In  this  direction,  the  work  recently  published  [42] 
showing  increased  mortality  of  TLR4  KO  mice  under  heat  stroke 
suggests  that  LPS  might  not  play  a  significant  pathogenic  role. 

In  summary,  the  present  work  provides  new  insights  into  the 
molecular  mechanisms  underlying  the  complex  etiology  of  HS  and 
defines  a  framework  that  supports  in  silico  exploration  of  cytokine 
signaling  pathways  in  response  to  HS.  This  type  of  modeling 
framework  not  only  aids  in  the  development  of  new  methods  that 
will  reduce  the  need  for  timely  and  costly  animal  experiments,  but 
also  increases  the  rapidity  and  accuracy  with  which  novel 
pharmacologic  intervention  and/or  treatments  are  identified  to 
treat  this  debilitating  illness.  However,  it  should  be  noted  that  this 
study  is  exploratory  and  a  larger  scale  study  is  needed  to  confirm 
the  results  in  the  future. 

Methods 

Ethics  Statement 

In  conducting  the  research  described  in  this  report,  the 
investigators  adhered  to  the  “Guide  for  Care  and  Use  of 
Laboratory  Animals”  as  prepared  by  the  Committee  on  Care 
and  Use  of  Laboratory  Animals  of  the  Institute  of  Laboratory 
Animal  Resources,  National  Research  Council.  The  protocol  was 
approved  by  the  Scientific  Review  Committee  and  the  Institu¬ 
tional  Animal  Care  and  Use  Committee  “US  Army  Research 
Institute  of  Environmental  Medicine  Institutional  Animal  Care 
and  Use  Committee”  (Permit  Number:  A09-02). 

Description  of  the  Data 

Details  of  the  HS  protocol  are  described  elsewhere  [43] .  In  this 
study,  conscious,  unrestrained  male  B6129F2  (wild-type;  WT)  and 
TNF-RI/R-II  knockout  (TNFR  KO)  mice  (originally  at  a  normal 
housing  temperature  of  25±0.2°C)  were  exposed  to  an  ambient 
temperature  (Ta)  of  39.5±0.2°C  in  an  incubator,  in  the  absence  of 
food  and  water,  until  a  Tcmax  of  42.4°C  was  attained.  Following 
removal  from  the  heat  at  Tc  max ,  food  and  water  were  provided  ad 
libitum  during  undisturbed  recovery  at  Ta  of  25±2°C. 

Prolonged  heat  exposure  induces  thermoregulatory  changes 
consisting  of  hyperthermia  in  response  to  direct  heat  exposure,  and 
a  biphasic  response  characterized  by  hypothermia  which  develops 
within  ~3  hours  of  recovery  [5,18].  The  initiation  of  the  SIRS  is 
thought  to  occur  within  the  time  frame  front  Tcnlax  to 
hypothermia,  during  which  endotoxin  is  thought  to  leak  across 
ischemic  damaged  gut  membranes  into  the  portal  circulation  [5,6]. 
To  identify  the  inflammatory  pathways  mediating  the  early  stages 
of  the  SIRS,  mice  were  assigned  to  one  of  the  following  groups  for 
blood  and  tissue  collection:  1)  baseline  (Tc<  36. 0cC;  immediately 
prior  to  heat  stress),  2)  TCymax  (Tc  =42.4°C),  3)  return  to  baseline 
(RTB;  first  Tc  value  <  36. 0°C  during  cooling),  or  4)  hypothermia 
depdi  (lowest  Tc  value  with  cooling  rate  of  O.OTC /min  during 
recovery).  Control  mice  were  tested  at  Ta  of  25°C  at  their  original 


cage  location  and  not  exposed  to  the  incubator  environment.  A 
HS  group  was  assigned  to  each  sampling  time  point,  except 
baseline,  which  was  represented  by  one  control  (nonheated)  group 
for  each  genotype.  Control  groups  were  included  at  each  time 
point  to  examine  circadian  influences  on  the  measured  variables  in 
the  absence  of  HS.  Group  sizes  for  each  sampling  time  point  were 
as  follows:  baseline:  WT,  n=  10;  TNFR  KO,  n  =  8;  Tcfnax\  WT 
control,  n  =  7;  TNFR  KO  control,  n  =:  6;  WT  heat,  n  =  7;  TNFR 
KO  heat,  n  =  6;  RTB:  WT  control,  n  =  8;  TNFR  KO  control, 
n  =  8;  WT  heat,  n  =  8;  TNFR  KO  heat,  n  =  9;  hypothermia:  WT 
control,  n  =  6;  TNFR  KO  control,  n  =  8;  WT  heat,  n  =  8;  TNFR 
KO  heat,  n  =  8.  Note  that  owing  to  the  need  for  sacrificing  mice  to 
obtain  biological  samples,  data  at  different  time  points  correspond 
to  different  mice,  which  introduced  additional  noise  due  to  the 
inherent  variability  among  animals.  Furthermore,  this  experimen¬ 
tal  design  precluded  an  analysis  of  survival  rates  in  this  study. 

Core  temperature  and  dehydration.  Tc  was  continuously 
monitored  at  1  -min  intervals  using  an  intraperitoneally  implanted 
battery-free  radiotelemetry  transmitter  (model  G2  Emitter,  Mini 
Mitter  Co.,  Inc.,  Bend,  OR)  in  conscious,  unrestrained  mice.  Mice 
were  implanted  with  a  transmitter  device  at  least  2  weeks  prior  to 
experimentation  to  ensure  full  recovery  prior  to  sample  collection. 
The  level  of  dehydration  was  estimated  by  the  difference  in  body 
weight  before  and  after  heat  exposure  determined  on  a  top-loading 
balance  (0.1  g)  and  corrected  for  transmitter  weight;  however,  we 
did  not  account  for  urine  of  feces  loss  [24],  Body  weight  at  the 
beginning  of  the  experiment  was  similar  between  genotypes  at 
~27  g. 

In  order  to  use  mean  values  for  the  genomic  and  proteomic  data 
obtained  from  different  animals,  an  average  of  the  Tc  profile  for 
each  of  the  mouse  strains  was  used  as  input  to  the  system  (see 
Figure  2).  The  sampling  points  for  data  collection  were  based  on 
mouse  Tc  responses;  therefore,  each  data  point  corresponding  to 
maximum  core  temperature  (TCj„wx),  RTB,  and  hypothermia 
depth  was  collected  at  a  different  time  point  of  recovery  with 
respect  to  the  onset  of  heat  exposure  (; i.e .,  mice  heat  and  cool  at 
different  rates  despite  similar  body  weights  prior  to  experimenta¬ 
tion).  For  this  reason,  averaging  along  the  time  axis  would  lead  to 
misleading  results  in  terms  of  temperature  values  (see  Supporting 
Information  SI);  thus  averaging  along  the  temperature  axis  was 
preferred.  To  this  aim,  we  proceeded  in  four  steps:  i)  filtered  the 
individual  curves,  ii)  divided  them  into  increasing  and  decreasing 
temperature  sections  (to  ensure  bijective  functions),  iii)  inversely 
mapped  by  interpolating  the  time  values  for  a  set  of  temperatures, 
and  iv)  computed  the  mean  with  respect  to  the  temperature  for  the 
interpolated  times.  Figure  2  shows  the  resultant  Tc  profiles  that 
were  used  as  an  input  to  the  mathematical  model.  Circadian  Tc, 
time  to  reach  Tcfnax ,  thermal  load,  cooling  time  to  baseline,  and 
time  to  hypothermia  were  analyzed  using  one-way  analysis  of 
variance  (AN OVA)  and  the  Holm-Sidak  method  for  post-hoc 
comparisons,  setting  the  significance  at  P<0.05. 

Previous  experiments  where  Tc  was  monitored  throughout  72 
hours  of  recovery  [43]  showed  that  after  hypothermia  (~6  hours) 
Tc  remains  at  normothermia  or  lower  until  ~24-32  hours  when  a 
fever-like  increase  is  observed.  The  model  is  based  on  the 
assumption  that  elevated  Tc  is  the  initiator  of  the  inflammatory 
response  by  means  of  DP,  LPS,  and  ROS;  therefore,  although 
temperature  measurements  after  hyperthermia  are  not  available 
for  this  experiment,  we  assumed  that  between  6  and  24  hours  no 
LPS  is  released  and  the  production  of  DP  and  ROS  is  not 
increased  from  its  basal  level. 

Liver  mRNA  accumulation.  The  threshold  cycle  (Ct)  for  a 
variety  of  genes  implicated  in  the  response  to  heat  stroke  was 
determined  in  liver  tissue  by  means  of  qPCR  (quantitative 
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polymerase  chain  reaction).  The  Ct  for  each  gene  was  defined  as 
the  PCR  cycle  at  which  the  emitted  fluorescence  signal  was  greater 
than  the  background  fluorescence  level  [44],  Amplification  was 
determined  to  be  detectable  if  Ct<  35.  Changes  in  mRNA 
accumulation  were  calculated  as  fold-change  relative  to  the 
controls  matched  to  the  same  sampling  point  using  the  2~AACt 
method  [45],  where 

AACf  =  AC'tff  ■—  ACtc  =  (Ctf/jarget  —  C't/fj/fc)  _ 

~(Ct c, target  ~  Ct C.Hk) 

being  CtHjarget,  CtHHK,  CtCjarget,  and  Ctc,HK  the  Ct  of  the 
target  and  housekeeping  (HK)  genes  under  HS  and  control 
conditions,  respectively.  qPCR  results  are  presented  with  asym¬ 
metric  standard  error  bars  in  the  figures  due  to  the  exponential 
relation  between  Ct  and  fold-change.  The  standard  deviation, 
propagating  the  error  of  both  control  and  experimental  groups,  is 
given  by 


°AA  ct  —  \J  a  Act  c  +  aA  C1H  (2) 

where  (Jaao  represents  the  standard  deviation  of  A  A  Ct;  o\Ct(, 
and  a\ctH  are  ^le  standard  deviation  of  the  difference  between 
target  and  HK  gene  threshold  cycles  for  mice  under  control  and 
heat  stress  conditions,  respectively. 

Liver  mRNA  accumulation  was  measured  for  HSP70,  NF-r'B 
family  members  (NF-kB1,  NF-rB2,  RelA,  RelB,  c-Rel,  and  IrcB), 
Jun  and  Fos  genes  (related  to  AP-1),  and  the  following  cytokines 
and  cytokine  receptors:  IL-la,  IL-l/I,  IL-1RI,  IL-1RII,  IL-6,  IL- 
6R,  gpl30,  IL-10,  TNF-a,  TNF-RI,  and  TNF-RII. 

Plasma  cytokines  data.  Plasma  level  of  cytokines  and 
soluble  receptors  was  determined  on  duplicate  samples  using  the 
FlowMetrix1M  System  (Luminex,  Austin,  TX),  which  permits  the 
simultaneous  quantitation  of  multiple  cytokines.  The  FlowMetrix 
System  is  a  “Multiplexed  Fluorescent  Bead-Based  Immunoassay”, 
with  the  kits  used  in  this  study  being  specific  for  mouse  cytokines. 
Sensitivities  of  the  cytokine  and  soluble  cytokine  receptor  assays 
were  ~2  and  ~20  pg/ml,  respectively.  The  plasma  cytokine  and 
soluble  receptor  data  used  to  calibrate  the  model  correspond  to  IL- 
lot,  IL-1  ft,  sIL-lRI,  sIL-lRII,  IL-6,  sIL-6R,  sgpl30,  IL-10,  TNF- 
a,  sTNF-RI,  and  sTNF-RII.  The  experimental  data  of  circulating 
cytokines,  including  control  and  HS  measurements,  were  analyzed 
to  detect  outliers  and  significant  changes  between  data  collected  at 
different  sampling  points.  The  error  bars  in  the  figures  represent 
the  standard  error  of  the  measures  (mean  +  SE).  The  original  data 
from  this  analysis  are  presented  elsewhere  [24], 

Modeling  and  Biological  Assumptions 

For  the  sake  of  parsimony,  we  used  first-order  kinetics  to  model 
activation  of  the  TFs  [46].  Therefore,  the  rate  of  activation  of  the 
TFs  is  proportional  to  the  concentration  of  the  activator  [Xj\  times 
the  concentration  of  the  inactive  TF,  denoted  by  [TF]: 

Rate  of  activation  =  /q [T;]  [7’F] .  (3) 


Moreover,  we  assumed  that  the  total  number  of  active,  [TFp], 
and  inactive  forms  of  TF  is  conserved: 


[TF0]  =  [TF]  +  [TFP],  (4) 


The  rate  of  change  of  [TFp]  is  given  by  the  balance  between  its 
activation  rate  by  all  the  possible  activators  (X(),  the  deactivation 
due  to  all  the  inhibitors  (Xj),  and  an  autonomous  deactivation  at  a 
rate  dpp: 


d[TFp\ 

dt 


^ activators 

E  k,[X,][TF] 

i=  1 


N inhibitors 

-  E  kj[Xj\[TFp\-dTF[TFp]. 
7  =  1 


(5) 


Since  the  signaling  networks  are  usually  very  complex  and 
information  is  lacking  from  our  study  regarding  the  intermediate 
states,  we  made  some  simplifications  through  lumped  states.  The 
activators  represent  the  ligand-receptor  complexes  or  other 
substances  as  DP  that  initiate  and/or  potentiate  the  signaling 
cascade. 

The  transcriptional  activation  was  modeled  using  Hill  functions 
and  mRNA  degradation  rate  was  assumed  to  be  linear: 


d[mRNAi] 

dt 


NTF,i 

E  kmRNAfJTFj 
7=1 


[TFj,p]n 


kk mRNAjJTFj  +  [^7,-P] 


-dlnRNAj[mRNAi\. 


(6) 


Mass  action  kinetics  was  used  for  describing  the  translation  of 
mRNA  into  proteins  (cytokines  and  cytokine  receptors)  and 
complex  formation: 

^  =  =  kprojj  mRNAj  [ill R N A  j]  dpro^[PvOtj] 

Nprotj 

-  E  kpwti_Proij[Proti][Protj] 

7=1  ’ 

Nprotj 

+  E  dpwti_Protj[ProtiJ,rotj\. 

7  =  1 

Note  that  kprotl^„p t/A,  are  not  pure  translation  rates  since  they 
also  account  for  the  transport  of  the  cytokines  outside  the  cell. 

The  soluble  receptors  involved  in  this  model  (sIL-lRI,  sIL-lRII, 
sIL-6R,  sgpl30,  sTNF-RI,  and  sTNF-RII)  are  produced  mainly  by 
proteolytic  cleavage  of  the  extracellular  domains  of  their  analogous 
membrane  bound  forms  [47 — 49] .  Therefore,  their  concentration 
was  modeled  as  follows: 


d[sProtj] 

^  =  ksProtj  [PYOti\  dgProti  \sPYOti\ 

^ sProti 

-  E  ksPwti_Protj[sProti\[ProtJ] 
7  =  1 

N sProti 

+  E  dsProti_Protj[sPrOti_PrOtj\ 
7  =  1 


(8) 
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and  the  corresponding  expression  for  the  analogous  cell-surface 
receptor  is  given  by  Eq.  (7)  plus  an  additional  term  accounting  for 
the  conversion  into  the  soluble  form  (  —  ksProt  [Proti]). 

Without  heat  the  model  is  assumed  to  be  at  steady  state 
(circadian  effects  were  not  statistically  significant): 

d[mRNA,]  _Q 
dt 


We  assumed  a  normalized  expression  rate  for  all  genes  and 
active  concentration  of  the  transcription  factors  to  be  equal  to  one 
( mRNAj  =  1  and  TFjP  =1).  Therefore,  we  computed  degradation 
rates  as: 


ntfj 
dmRNAj  =  Y. 


7=1 


kmRNA, 


kkn  4-  1 

KKmRNAj_TFj  ^  1 


(10) 


and  significantly  decreased  the  number  of  parameters  to  be 
estimated  (from  217  to  130). 

LPS.  Different  types  of  environmental  stressors,  including 
heat,  can  cause  increased  intestinal  permeability  that  facilitates 
endotoxin  leakage  across  ischemic-damaged  gut  epithelial  mem¬ 
branes  and  induces  local  and/ or  systemic  inflammatory  reactions 
[50].  The  presence  of  circulating  endotoxin  in  some  HS  patients  is 
thought  to  be  due  to  leakage  following  gastrointestinal  barrier 
disruption,  but  may  also  be  related  to  liver  damage  [5 1] .  Since  the 
liver  is  one  of  the  major  sites  of  endotoxin  clearance,  damage  to 
this  organ  may  result  in  increased  susceptibility  of  mammals  to  the 
heat-induced  SIRS  mediated  by  endotoxin. 

Su  et  al.  [52]  performed  in  vivo  analyses  of  physiologically 
relevant  barrier  dysfunction  and  determined  mouse  jejunal 
permeability  by  measuring  the  paracellular  bovine  serum  albumin 
(BSA)  flux  under  standard  conditions  as  0.18  [iBSA/(h  *  cm). 
Lambert  et  al.  [12]  and  Dokadny  et  al.  [53]  showed  significant 
time-  and  temperature-dependent  increases  in  gastrointestinal 
permeability  following  both  modest  and  severe  heat  insults. 
Although  several  cytokines  and  other  factors  are  known  to  regulate 
the  permeability  of  the  gut  epithelial  barrier  [54] ,  in  this  work  we 
only  considered  permeability  influenced  by  temperature,  assuming 
an  exponential  increase  front  the  basal  value  (0.18 
/ tBSA/(h  *  cm))  when  Tc  is  above  41°C.  Increased  gut  epithelial 
barrier  permeability  at  this  core  temperature  resulted  in  the  release 
of  LPS  which  was  modeled  using  mass  action  kinetics  with  the 
following  equation: 


=  kLPS  ( [LPS^]  -  [LPS\ )  -  dLPS  [LPS]  (11) 


cytokines,  chemokines,  and  adhesion  molecules  [57].  Moreover, 
oxidative  stress  impairs  the  heat  stress  response  and  delays 
unfolded  protein  recovery  and  function,  which  may  compromise 
protective  functions  of  some  protein  pathways  during  heat  stress 

[58] . 

Zhang  et  al.  [36]  observed  a  small,  transient  increase  in  hepatic 
oxidative  damage  in  young  rats  that  experienced  minimal  liver 
damage  under  heat  stress;  however,  this  response  was  accompa¬ 
nied  by  a  sharp  and  significant  elevation  of  AP-1  DNA  binding 
activity  suggesting  that  this  TF  is  mediating  inflammatory  changes 
in  the  liver  during  recovery.  In  contrast,  NF-;cB  DNA  binding 
activity  did  not  change  following  heat  stress  indicating  that  NF-jcB 
might  not  be  a  main  regulator  of  cytokine  gene  expression  under 
these  conditions.  It  should  be  noted  that  the  results  of  this  study 
cannot  be  examined  solely  in  the  context  of  heat  stress  because  the 
research  design  did  not  control  for  dehydration  that  was 
experienced  during  heat  exposure,  which  is  an  adverse  physiolog¬ 
ical  response  that  is  known  to  increase  oxidative  and  cellular  stress 

[59] .  In  order  to  encompass  these  experimental  studies,  ROS  are 
generated  in  our  model  following  zero-order  kinetics  with  varying 
kRos  and  degraded  following  first-order  kinetics: 

=kRos~dRos[BOS]  (12) 

where  kRos  is  a  function  of  the  temperature  and  the  dehydration 
level.  ROS  activates  AP-1  transcription  following  equation  (5)  and 
inhibits  HSF  as  described  below  in  equation  (18). 

Denatured  proteins,  HSPs,  and  HSF-1.  HSPs  are  impor¬ 
tant  modulators  of  both  anti-inflammatory  and  pro-inflammatory 
responses.  Asea  et  al.  [60]  introduced  the  term  chaperokine,  to 
describe  the  dual  role  of  most  HSPs  as  chaperones  with  an 
intracellular  cytoprotective/antiapoptotic  function  and  cytokines 
with  an  extracellular  immunogenic  function.  In  non-stressed 
conditions,  HSPs  function  as  molecular  chaperones  by  maintain¬ 
ing  protein  conformation  and  facilitating  transport  throughout  the 
cell’s  various  compartments.  Under  environmental  stress  condi¬ 
tions,  the  structural  integrity  of  cellular  proteins  is  compromised 
and  binding  of  HSPs  to  these  damaged  proteins  prevents 
aggregation  and  supports  refolding. 

Hyperthermic  temperatures  cause  accumulation  of  denatured 
proteins  and  exposure  of  their  hydrophobic  domains  which 
stimulates  the  expression  of  HSP  genes  (especially  HSP70,  the 
most  highly  heat-inducible  member)  through  activation  of  HSF-1 
[61].  The  onset  temperature  for  this  response  is  usually  ~40t>C, 
but  some  transitions  may  extend  as  low  as  37  — 38cC  [62]; 
therefore,  we  utilized  the  equation  in  [63]  to  calculate  the 
fractional  protein  denaturation,  Vden,  as  a  function  of  the 
temperature  in  the  range  of  37  — 45°C: 


where  kPPs  represents  the  gut  epithelial  barrier  permeability  and 
increases  with  the  core  temperature,  [ LPSguts }  is  the  initial 
concentration  of  LPS  in  the  gastrointestinal  tract  of  the  mouse  and 
[LPS]  is  the  concentration  of  LPS  in  the  circulation.  LPS  acts  as 
the  prototypical  endotoxin  because  it  binds  the  Toll-like  receptor  4 
(TLR4)  to  induce  a  signal  transduction  cascade  that  ultimately 
triggers  essential  signaling  modules  resulting  in  activation  of  the 
transcription  factors  NF-kB  and  AP-1  [55]  following  equation  (5). 

Oxidative  Stress.  Heat  stress  increased  reactive  oxygen 
species  (ROS)  generation  and  free  radical-mediated  splanchnic 
injury  of  young  rats  [56].  Furthermore,  reactive  oxygen  and 
nitrogen  species  production  is  thought  to  be  involved  in  regulation 
of  redox-sensitive  transcription  factors,  such  as  AP-1  and  NF-kB, 
that  mediate  the  expression  of  inflammatory  mediators  such  as 


Hrfe„=(^l_-04jx0.03xL4r-37.  (13) 

The  amount  of  protein  that  denatures  per  unit  of  time  is 
obtained  by  multiplying  Vden  by  the  amount  of  native  proteins  (P): 

Qden=Vde„[P}.  (14) 

The  native  proteins  are  obtained  by  subtracting  the  denatured 
proteins  (DP)  front  the  total  amount  of  native  proteins  susceptible 
to  denaturation  in  the  range  of  37-45°C  (P0),  approximated  as 
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10%  of  the  total  amount  of  proteins  in  the  cell: 

[P0]  =  [P]  +  [DP]  +  [HSPMP] .  (15) 

Consistent  with  [63],  the  amount  of  denatured  proteins  which 
complex  with  HSP70  is  determined  by  mass  balance  in  the  present 
model.  In  contrast  to  the  other  mass  balances  in  the  model  in 
which  there  is  an  interaction  between  two  substances  to  form  a 
complex,  in  this  mass  balance  free  HSP70  interacts  with  both  the 
denatured  and  the  native  proteins:  free  HSP70  ( HSP )  binds  to  free 
DP  (DP),  reducing  the  amount  of  free  DP  in  the  cell,  while  the 
released  (renatured)  proteins  are  added  to  the  pool  of  native 
proteins  in  the  cell.  The  complexation  is  described  by  the  following 
relation: 

d\IISP  DP]  _ 
dt 

kHSP_DP[HSP\[DP ]  +kHSF  hsp  dp[HSFMSP][DP]  (16) 

-dHSp_Dp[HSP-DP]. 


The  first  term  describes  the  forming  of  the  HSP_DP  complexes 
while  the  last  term  describes  the  dissociation  of  the  complexes  into 
renatured  proteins.  The  amount  of  free  denatured  proteins,  DP,  is 
obtained  from  the  amount  of  denaturing  proteins  per  unit  of  time, 
Qden,  minus  the  amount  of  DP’s  complexed  with  HSP70: 


d[DP] 

dt 


Qden  ~  kHSP_Dp[HSP]  [DP] 


—  knsF_HSP_DP  [HSFMSP]  [DP] . 


(17) 


Before  heat-induced  activation,  HSF-1  exists  as  a  monomer 
localized  to  the  cytoplasm  with  much  of  it  forming  a  heterodimer 
complex  with  HSP70.  HSP70  preferentially  binds  to  denatured 
proteins;  therefore,  it  has  been  postulated  that  activation  of  HSF-1 
may  occur  as  a  result  of  competitive  release  of  this  TF  from  the 
HSF1-HSP70  complex  when  the  concentration  of  denatured 
cytoplasmic  proteins  increases  as  a  result  of  heat  shock  [26,64],  As 
a  result,  the  equation  describing  HSF-1  dynamics  reads: 

=  kjjsF-HSF-DP  [HSF  HSP]  [DP]  -  kHSF_HSP  [HSF]  [HSP] 

(18) 

-  kHSF_ROs  [HSF]  [ROS]  +  dHSF_HSP  [HSF  MSP] 

which  is  equivalent  to  eq.  (5)  being  [HSF MSP]  the  inactive  state 
of  HSF,  [DP]  the  activator,  and  [ROS]  the  inhibitor. 

HSP  mRNA  (mHSP)  transcription  was  modeled  using  eq.  (6) 
where  the  only  TF  is  HSF.  The  Hill  coefficient  for  HSF-1  was 
taken  from  [65]  as  1.5.  The  intermediate  value  for  the  apparent 
Hill  coefficient  of  between  1  and  2  suggests  that  one  HSF  trimer 
may  bind  stably,  and  a  second  may  bind  only  weakly  or  partially. 
mHSP  is  then  translated  into  HSPs  through  eq.  (7),  which  exert  an 
anti-inflammatory  effect  by  inhibiting  translocation  of  NF-rB  to 
the  nucleus  and  preventing  expression  of  inflammatory  mediators 
in  a  dose-  and  time-dependent  manner  [35].  Moreover,  the  crucial 
indirect  role  of  HSPs  in  maintaining  gut  epithelial  barrier  integrity 
suggests  an  important  anti-inflammatory  effect  by  attenuating 
endotoxin  leakage  into  the  circulation,  which  should  mitigate  the 
heat-induced  SIRS  [66]. 


NF-rB.  NF-rB  is  a  ubiquitous  TF  of  particular  importance  in 
immune  and  inflammatory  responses.  The  larger  NF-rB  family  is 
composed  of  two  subfamilies:  the  NF-rB  subfamily  that  includes 
NF-kB1  (p50/pl05)  and  NF-rB2  (p52/pl00)  proteins,  and  the 
Rel  subfamily  that  includes  RelA  (p65),  RelB,  and  c-Rel.  All  these 
structurally-related  proteins  can  form  homodimers  or  heterodi- 
mers  (except  for  RelB,  which  only  forms  heterodimers).  NF-rB 
belongs  to  the  class  of  “rapid-acting”  primary  TFs,  which  are 
present  in  cells  in  an  inactive  state  and  do  not  require  protein 
synthesis  to  be  activated.  Rel/NF-tcB  transcription  complexes  are 
present  in  a  latent,  inactive  state  in  the  cytoplasm  due  to  binding  to 
the  inhibitor  IrB.  A  variety  of  stimuli,  namely  LPS  through 
binding  to  TLR4,  rapidly  activate  Rel/NF-R'B  transcription 
complexes  by  releasing  the  TFs  from  inhibitor  binding,  which 
allows  translocation  to  the  nucleus  for  binding  to  rB  sites  that 
regulate  the  expression  of  many  genes  [67] .  The  pro-inflammatory 
cytokines  IL-la,  IL-l/f,  and  TNF-ot  experience  reciprocal  activa¬ 
tion  with  NF-rB  [68].  NF-rB  also  regulates  the  expression  of 
additional  genes  encoding  proteins  involved  in  the  heat  shock 
response,  including  IL-IRa,  IL-1RII,  IL-6,  IL-6R,  gp-130,  IL-10, 
IL-10R,  and  TNF-RII  [27,68-70], 

The  development  of  mathematical  models  of  NF-jcB  signaling, 
tightly  linked  to  experimental  results,  has  been  instrumental  in 
unraveling  the  forms  of  regulation  in  NF-rB  signaling  and  their 
underlying  molecular  mechanisms  [71,72].  Unfortunately,  due  to 
the  lack  of  data  on  nuclear  translocation,  we  were  not  able  to 
incorporate  details  on  NF-rB  activation  into  our  model;  dius,  we 
assumed  first  order  kinetics  (Eq.  5)  treating  the  aforementioned  TF 
as  a  lumped  state.  The  apparent  Hill  coefficient  for  the 
transcriptional  activation  of  downstream  genes  was  taken  from 
the  literature  and  assumed  to  have  a  value  of  2  [73,74], 

AP-1.  AP-1  is  another  essential  TF  that  regulates  inflamma¬ 
tory  and  immune  genes.  AP-1  is  a  group  of  structurally  and 
functionally  related  members  of  the  Jun  and  Fos  protein  families. 
Jun  proteins  exist  as  homo-  and  heterodimers  whereas  Fos 
proteins,  which  cannot  homodimerize,  form  stable  heterodimers 
with  Jun  proteins,  thereby  enhancing  their  DNA-binding  activity. 
AP-1  activity  is  regulated  by  a  broad  range  of  physiological  and 
pathological  stimuli,  including  cytokines,  growth  factors,  stress 
signals,  infectious  agents,  and  oncogenic  stimuli.  Regulation  of  net 
AP- 1  activity  can  be  achieved  through  changes  in  transcription  of 
genes  encoding  AP-1  subunits,  stability  of  their  mRNAs, 
posttranslational  processing  and  turnover  of  pre-existing  or  newly 
synthesized  AP- 1  subunits,  and  specific  interactions  between  AP- 1 
proteins  and  other  TFs  and  cofactors  [75], 

Developing  a  detailed  model  for  AP-1  activation  is  out  of  the 
scope  of  this  work;  therefore,  we  assumed  first  order  kinetics  (Eq. 
5)  considering  AP-1  as  a  lumped  state  activated  by  LPS  through 
TLR4,  IL-la,  11,-1  /i,  IL-6,  and  ROS.  Among  the  genes  involved 
in  the  cytokine  network  with  HS,  nine  are  known  to  be  AP-1  target 
genes:  IL-la,  IL-lft  IL-IRa,  IL-1RI,  IL-1RII,  IL-6,  IL-10,  TNF- 
RI,  TNF-RII  [26,76,77].  Transcriptional  activation  of  these  genes 
was  modeled  using  Eq.  (6)  with  a  Hill  coefficient  of  2  (AP- 1  Hill 
coefficient  was  found  to  vary  between  1.6  and  2.6  [78]). 

STAT-3.  STAT-3  is  a  TF  with  fundamental  importance  for 
cytokine-mediated  induction  of  acute-phase  response  genes  and  is 
a  key  regulator  of  gene  expression  in  response  to  IL-10  and 
glycoprotein  130  (gp  1 30)  family  cytokine  signaling  (e.g.,  IL-6) 
[79,80].  IL-10  acts  as  a  more  potent  anti-inflammatory  cytokine 
than  IL-6,  although  both  cytokines  activate  STAT-3.  Cytokine 
binding  to  cell  surface  receptors  induces  receptor-associated  Janus 
tyrosine  kinase  1  JAK1)  activation  leading  to  phosphorylation  of  a 
single  tyrosine  residue  in  the  STAT-3  molecule.  Phosphorylation 
results  in  STAT-3  dimerization  and  nuclear  entry  for  binding  to 
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specific  DNA  sequences  in  the  promoter  regions  of  target  genes 
[80]. 

STAT-3  is  thought  to  upregulate  the  following  proteins  involved 
in  the  cytokine  network  during  heat  stroke:  IL-1R1,  IL-6R,  gpl30, 
IL-10,  and  suppressor  of  cytokine  signaling  3  (SOCS-3)  [70,80 — 

82] .  The  initial  activation  of  STAT-3  by  IL-10  or  IL-6  (Eq.  5) 
precedes  SOCS-3  gene  expression  (Eq.  6)  with  subsequent 
inhibitory  effects  of  SOCS-3  on  gpl30/IL-6  signaling  padiways 

[83]. 

IL-1  family.  IL- 1  is  an  important  mediator  of  inflammation 
and  tissue  damage  in  multiple  organs.  The  IL-1  family  consists  of 
two  agonists  (IL-1  a  and  IL-1 /I),  two  receptors  (biologically  active 
IL-1RI  and  inert  IL-1RII),  and  a  specific  receptor  antagonist  (IL- 
lRa).  IL-1  signals  through  a  single  receptor  (IL-1RI)  that  binds  IL- 
la  and  IL-1  fi  with  equal  affinity,  while  IL-1RII  preferentially 
binds  IL-1  p.  IL- 1  Ra  prevents  the  association  of  the  IL- 1  ligands 
with  the  IL-IRs;  thus,  the  balance  between  IL-1  and  IL-IRa  in 
local  tissues  plays  an  important  role  in  susceptibility  and  severity  of 
many  disease  conditions  [27,84],  IL-1  signaling  activates  NF-kB 
and  AP-1  to  induce  the  expression  of  many  genes  [85]. 

Both  IL-1RI  and  IL-1RII  have  soluble  isoforms  (sIL-lRI  and 
sIL-lRII)  that  are  generated  by  proteolytic  cleavage  of  their 
extracellular  domains  modeled  by  Eq.  (8).  These  sIL-IRs  bind  and 
inhibit  IL- 1  a  and  IL- 1  /?  signaling.  Because  IL- 1  Ra  binds  avidly  to 
both  the  membrane-bound  and  soluble  forms  of  IL-1  RI,  sIL-lRI 
is  likely  the  principal  antagonist  of  IL-IRa.  Therefore,  although 
sIL-lRII  antagonizes  the  action  of  IL-1/?,  sIL-lRI  indirectly  may 
serve  to  facilitate  the  activity  of  IL-la  and  IL- 1  /?  by  binding  to  IL- 
IRa  [86], 

IL-6  family.  Circulating  IL-6  shows  the  highest  correlation 
with  mortality  and  neurologic  symptoms  in  HS  patients,  suggest¬ 
ing  this  cytokine  may  be  an  important  therapeutic  target  for 
prevention/treatment  strategies.  However,  IL-6  KO  mice  showed 
higher  mortality  rates  than  WT  mice  indicating  that  IL-6  also  has 
protective,  anti-inflammatory  actions  that  are  critical  for  survival 
[18].  IL-6  is  an  important  mediator  of  the  acute-phase  response  to 
injury  and  infection  and  induces  its  cellular  actions  through  two 
signaling  pathways  that  have  opposing  actions  [87].  In  the 
“classic”  (anti-inflammatory)  pathway,  IL-6  first  binds  to  its  non¬ 
signaling  membrane-bound  IL-6R  (also  called  a-receptor  subunit, 
IL-6Ra)  followed  by  recruitment  of  the  signaling  transducing 
receptor  subunit  gpl30  to  the  complex,  resulting  in  activation  of 
anti-inflammatory  cascades  [88,89].  Contrasting  this,  the  IL-6 
trans-signaling,  pro-inflammatory  pathway  is  activated  when  IL-6 
binds  to  the  soluble  isoform  sIL-6R  and  forms  complexes  that 
intercalate  into  the  membranes  of  cells  that  contain  gpl30,  but 
normally  do  not  respond  to  the  cytokine  [48] .  The  trans-signaling 
pathway  is  modulated  by  sgpl30,  a  circulating  cleavage  product  of 
the  membrane-bound  receptor  subunit  [90].  In  our  model,  the 
dual  opposing  actions  of  IL-6  were  captured  by  assuming  that  the 
“classic”  pathway  activates  the  anti-inflammatory  transcription 
factor  STAT-3  and  the  trans-signaling  pathway  activates  AP-1 .  IL- 
6  strongly  induces  SOCS-3  protein  through  STAT-3  and  in  turn, 
IL-6  signaling  is  selectively  inhibited  owing  to  the  binding  of 
SOCS-3  to  the  IL-6R  subunit  gpl30  [91], 

IL-10  family.  IL-10  is  a  potent  anti-inflammatory  cytokine 
that  modulates  the  expression  of  several  cytokines,  including  IL- 1 , 
IL-6  and  TNF-a  [76].  A  major  function  of  IL-10  is  to  control  and 
reduce  excessive  immune  responses  during  infection  and  autoim¬ 
munity,  mainly  by  inhibiting  the  production  of  pro-inflammatory 
cytokines  in  macrophages  and  other  cell  types  [92,93].  IL-10  can 
induce  the  expression  of  SOCS-3,  suggesting  that  the  capacity  of 
IL- 1 0  to  inhibit  the  expression  of  LPS-inducible  pro-inflammatory 
genes  may  depend  on  SOCS-3  [91],  In  this  model,  IL-10- 
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induction  of  SOCS-3  was  mediated  by  the  transcription  factor 
STAT-3. 

TNF  family.  TNF-a  is  a  pro-inflammatory  cytokine  with 
important  actions  in  immunity  and  inflammation,  including  the 
control  of  cell  proliferation,  differentiation,  and  apoptosis.  Binding 
of  TNF-o!  to  its  two  receptors,  TNF-RI  and  TNF-RII,  results  in 
recruitment  of  signal  transducers  that  activate  at  least  three  distinct 
effectors.  Through  complicated  signaling  cascades  and  networks, 
these  effectors  activate  caspases  as  well  as  AP-1  and  NF-kB 
[94,95].  Generation  of  soluble  TNF-RI  and  TNF-RII,  by 
proteolytic  cleavage,  is  also  considered  a  highly  regulated  process. 
These  circulating  soluble  receptors  modify  ligand  actions  by 
stabilizing  TNF-a  protein  structure,  decreasing  membrane  recep¬ 
tor  number,  or  specifically  inhibiting  ligand-receptor  binding  [86] . 

Model  Calibration 

Model  calibration,  or  parameter  estimation,  is  a  key  step  in  the 
development  of  reliable  dynamic  models.  Given  a  model  structure 
and  a  set  of  experimental  data,  the  objective  of  parameter 
estimation  is  to  calibrate  the  model  to  reproduce  the  experimental 
results  in  the  best  possible  way.  It  is  usually  formulated  as  the 
optimization  of  a  scalar  cost  function,  J(p),  which  measures  the 
goodness  of  the  fit  with  respect  to  the  model  parameters  peU.NF . 
This  function  consists  of  a  weighted  distance  measure  between  the 
experimental  values  corresponding  to  the  measured  variables, 
represented  by  the  vector  y,  and  the  predicted  values  for  those 
variables,  represented  by  the  vector  y.  Several  estimator  functions 
have  been  suggested  as  metrics,  where  the  weighted  least-squares 
estimator  is  the  most  common  [96]: 

NENV/NMij 

•4(p)  wyk \yi]k -yijk( p)  •  (19) 

1=1 j=\  k=\ 

Here,  NE  is  the  number  of  experiments,  NVf  number  of 
measured  outputs  in  experiment  i,  NMy  number  of  measurements 
of  output  j  during  experiment  i,  yyk  model  predicted  value  k  of 
output  j  in  experiment  i,  yyk  measurement  k  of  output  j  in 
experiment  i,  and  Wyk  the  weight  of  measurement  k  of  output  j  in 
experiment  i. 

Special  attention  must  be  paid  to  the  selection  of  the  weights 
since  the  optimal  value  of  p  will  depend  on  them.  When  a  good 
approximation  for  the  standard  deviation  of  the  data  is  available,  a 
good  choice  for  the  weights  is  Wyk  =  1  / ojjk  where  a yk  is  the 
standard  deviation  of  measurement  k  of  output  j  in  experiment  i. 
In  this  case,  minimizing  J/s  is  equivalent  to  minimizing  the 
Maximum  Likelihood  Estimator  introduced  by  Fisher  in  1912 
[97],  which  maximizes  the  probability  of  the  observed  event. 
However,  in  the  present  study,  preliminary  fitting  indicated  that 
for  some  of  the  species  there  is  no  combination  of  parameters  that 
accurately  fits  the  liver  and  plasma  data  simultaneously.  In  order 
to  test  the  hypothesis  that  the  liver  is  one  of  the  major  sources  of 
circulating  cytokines,  we  prioritized  the  fit  of  the  liver  mRNA 
accumulation  by  increasing  the  weights  of  the  qPCR  data  with 
respect  to  those  of  the  soluble  cytokines. 

Due  to  the  nonlinear  nature  of  the  model  considered  here,  the 
resulting  optimization  problem  is  multimodal  (non-convex). 
Therefore,  traditional  gradient  based  methods,  like  Levenberg- 
Marquardt  or  Gauss-Newton,  may  fail  to  identify  the  global 
solution  and  may  converge  to  a  local  minimum  when  a  better 
solution  exists  just  a  small  distance  away.  Moreover,  in  the 
presence  of  a  bad  fit,  there  is  no  way  of  knowing  if  it  is  due  to  a 
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wrong  model  formulation,  or  if  it  is  simply  a  consequence  of  local 
convergence.  Thus,  there  is  a  distinct  need  for  using  global 
optimization  methods  which  provide  more  guarantee  of  converg¬ 
ing  to  die  globally  optimal  solution.  In  this  work,  the  model  was 
fitted  to  the  available  experimental  data  using  SSm,  a  global 
optimization  metaheuristic  based  on  Scatter  Search  developed  for 
parameter  estimation  in  nonlinear  dynamic  biochemical  systems 
[98]  and  available  in  the  toolbox  SensSB  [22]. 

Sensitivity  and  Correlation  Analysis 

A  practical  identifiability  analysis  aims  to  determine  whedier, 
given  a  model  structure,  the  parameters  of  a  model  could  be 
uniquely  identified  front  the  available  (limited  and  noisy)  data  [99] . 
There  are  two  main  aspects  that  influence  model  identifiability: 
the  sensitivity  of  the  parameters  and  the  correlation  among  them. 

Sensitivity  analysis  indicates  which  parameters  are  the  most 
important  and  therefore  would  have  the  greatest  impact  on  the 
predictions  of  the  model.  To  analyze  how  the  model  variables 
change  around  the  best  parameter  set  obtained,  we  computed 
local  sensitivity  coefficients  that  are  the  partial  derivatives  of  the 
model  state  variables  to  the  model  parameters  evaluated  at  the 
optimal  point.  To  make  these  measures  comparable  for  param¬ 
eters  and  states  of  different  order  of  magnitude,  relative  measures 
were  used  where  the  sensitivity  function  is  normalized  by  the  value 
of  the  parameter  and  the  state: 


s=Pq( 

yj  '■^Pq)  y=y(t,p),p  =  p 


(20) 


The  sensitivity  of  all  the  measured  states  with  respect  to  one 
parameter  can  be  summarized  as: 


^ msqr _ 


(21) 


Thus,  a  large  value  of  S"sqr  would  mean  that  a  change  in  the 
parameter  pq  has  an  important  effect  on  the  model  outcome.  This 
makes  the  parameter  identifiable  with  the  data  available  if  all  the 
other  parameters  are  fixed  and  the  larger  the  sensitivity,  the  more 
accurately  a  single  parameter  can  be  identified.  Therefore,  values 
of  critical  parameters  can  be  refined  while  parameters  having  little 
effect  can  be  simplified  or  even  ignored. 

Although  necessary,  high  parameter  sensitivity  is  not  sufficient 
to  ensure  the  identifiability  of  the  model.  In  the  case  of  several 
parameters,  the  sensitivity  functions  of  the  parameters  have  to  be 
linearly  independent.  The  degree  of  linear  dependence  among  the 
sensitivity  functions  can  be  measured  by  means  of  a  correlation 
analysis  based  on  the  Fisher  Information  Matrix  (FIM)  as 
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described  in  [100].  Correlations  among  parameters  close  to  +1 
or  —  1  mean  that  the  parameters  are  not  individually  identifiable 
because  a  change  in  one  parameter  can  be  compensated  by 
changes  in  the  other  parameters.  In  that  case,  an  infinite  number 
of  parameter  sets  fitting  the  experimental  data  with  the  same 
accuracy  would  exist,  thus  making  the  confidence  intervals  very 
large.  For  this  reason,  the  model  should  be  reduced  by  fixing  some 
of  the  parameters  to  their  nominal  values  or  by  properly  grouping 
some  sets. 
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Supporting  Information  S2  Schematic  diagram  of  the  cellular 
network  of  interactions  amongst  HSP70,  TLR4,  IL-1,  IL-6,  IL-10, 
and  TNF  families  induced  by  heat  stroke. 
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Supporting  Information  S3  Model  equations  and  SensSB  files 
needed  to  reproduce  the  results. 

(XLSX) 

Supporting  Information  S4  Table  with  the  list  of  parameters 
and  their  values. 

(ZIP) 

Supporting  Information  S5  Matlab  figures  with  the  simulated 
results  for  all  the  variables  for  WT,  TNFR  KO,  IL-10  KO,  high 
temperature,  and  LPS  injection  experiments. 

(ZIP) 
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